% plot December 2011 minus November 2011 of PhiBot, Eta, ThetaBot, Theta
clear all, close all
cd ~dmenemen/llc_4320/regions/Kuroshio
nx=3662; ny=392; nz=90;

% PhiBot
PhiBotNov=zeros(nx,ny); N=0;
for t=datenum(2011,11,1):(1/24):datenum(2011,11,30,23,0,0)
    fnm=['PhiBot/PhiBot_3662x392.' datestr(t,30)];
    tmp=readbin(fnm,[nx ny]);
    if mean(tmp)>500
        PhiBotNov=PhiBotNov+tmp;
        N=N+1;
    else
        disp(minmax(tmp))
    end
end
PhiBotNov=PhiBotNov/N;
PhiBotDec=zeros(nx,ny); N=0;
for t=datenum(2011,12,1):(1/24):datenum(2011,12,31,23,0,0)
    fnm=['PhiBot/PhiBot_3662x392.' datestr(t,30)];
    tmp=readbin(fnm,[nx ny]);
    if mean(tmp)>500
        PhiBotDec=PhiBotDec+tmp;
        N=N+1;
    else
        disp(minmax(tmp))
    end
end
PhiBotDec=PhiBotDec/N;

clf
subplot(311)
mypcolor(PhiBotNov');
thincolorbar
subplot(312)
mypcolor(PhiBotDec');
thincolorbar
subplot(313)
mypcolor(PhiBotDec'-PhiBotNov');
thincolorbar



% Eta
EtaNov=zeros(nx,ny); N=0;
for t=datenum(2011,11,1):(1/24):datenum(2011,11,30,23,0,0)
    fnm=['Eta/Eta_3662x392.' datestr(t,30)];
    EtaNov=EtaNov+readbin(fnm,[nx ny]);
    N=N+1; mydisp(N)
end
EtaNov=EtaNov/N;
EtaDec=zeros(nx,ny); N=0;
for t=datenum(2011,12,1):(1/24):datenum(2011,12,31,23,0,0)
    fnm=['Eta/Eta_3662x392.' datestr(t,30)];
    EtaDec=EtaDec+readbin(fnm,[nx ny]);
    N=N+1; mydisp(N)
end
EtaDec=EtaDec/N;

% ThetaBot
ThetaBotNov=zeros(nx,ny); N=0;
for t=datenum(2011,11,1):(1/24):datenum(2011,11,30,23,0,0)
    fnm=['ThetaBot/ThetaBot_3662x392.' datestr(t,30)];
    ThetaBotNov=ThetaBotNov+readbin(fnm,[nx ny]);
    N=N+1; mydisp(N)
end
ThetaBotNov=ThetaBotNov/N;
ThetaBotDec=zeros(nx,ny); N=0;
for t=datenum(2011,12,1):(1/24):datenum(2011,12,31,23,0,0)
    fnm=['ThetaBot/ThetaBot_3662x392.' datestr(t,30)];
    ThetaBotDec=ThetaBotDec+readbin(fnm,[nx ny]);
    N=N+1; mydisp(N)
end
ThetaBotDec=ThetaBotDec/N;

% Theta
ThetaNov=zeros(nx,ny,nz); N=0;
for t=datenum(2011,11,1):(1/24):datenum(2011,11,30,23,0,0)
    fnm=['Theta/Theta_3662x392x90.' datestr(t,30)];
    ThetaNov=ThetaNov+readbin(fnm,[nx ny nz]);
    N=N+1; mydisp(N)
end
ThetaNov=ThetaNov/N;
ThetaDec=zeros(nx,ny,nz); N=0;
for t=datenum(2011,12,1):(1/24):datenum(2011,12,31,23,0,0)
    fnm=['Theta/Theta_3662x392x90.' datestr(t,30)];
    ThetaDec=ThetaDec+readbin(fnm,[nx ny nz]);
    N=N+1; mydisp(N)
end
ThetaDec=ThetaDec/N;

% depth average Theta
region_name='Kuroshio';
pin=['~dmenemen/llc_4320/regions/' region_name '/'];
hFacC=readbin([pin 'grid/hFacC_3662x392x90'],[3662,392,90]);
load ../WPac/grid/thk90
for k=1:90
    hFacC(:,:,k)=hFacC(:,:,k)*thk90(k);
end
ThetaNovInt=sum(ThetaNov.*hFacC,3)./sum(hFacC,3);
ThetaDecInt=sum(ThetaDec.*hFacC,3)./sum(hFacC,3);

EtaDec      = single(EtaDec);
EtaNov      = single(EtaNov);
PhiBotDec   = single(PhiBotDec);
PhiBotNov   = single(PhiBotNov);
ThetaBotDec = single(ThetaBotDec);
ThetaBotNov = single(ThetaBotNov);
ThetaDec    = single(ThetaDec);
ThetaDecInt = single(ThetaDecInt);
ThetaNov    = single(ThetaNov);
ThetaNovInt = single(ThetaNovInt);

clear N a* d* f* n* t* h* k p* r*
save figs/NovDec

%%%%%%%%%%%%%%%

clear all, close all
cd ~dmenemen/llc/llc_4320/regions/Kuroshio/figs
load NovDec
orient landscape
wysiwyg

clf
subplot(411)
mypcolor((PhiBotDec'-PhiBotNov')/9.81*100);
set(gca,'xticklabel',[],'yticklabel',[])
caxis([-1 1]*6)
thincolorbar
title('November to December 2011 bottom pressure change (cm)')
subplot(412)
mypcolor((EtaDec'-EtaNov')*100);
set(gca,'xticklabel',[],'yticklabel',[])
caxis([-1 1]*25)
thincolorbar
title('November to December 2011 sea surface height change (cm)')
subplot(413)
mypcolor((ThetaBotDec'-ThetaBotNov'));
set(gca,'xticklabel',[],'yticklabel',[])
caxis([-1 1]/40)
thincolorbar
title('November to December 2011 bottom temperature change (^oC)')
subplot(414)
mypcolor((ThetaDecInt'-ThetaNovInt'));
set(gca,'xticklabel',[],'yticklabel',[])
caxis([-1 1]/4)
thincolorbar
title('November to December 2011 depth-avereage temperature change (^oC)')
print -dpng NovDec

%%%%%%%%%%%%%%%
% plot time series
clear all, close all
cd ~dmenemen/llc_4320/regions/Kuroshio
nx=3662; ny=392; nz=90;
dt=datenum(2011,11,1):(1/24):datenum(2011,12,31,23,0,0);
nt=length(dt);

% PhiBot
PhiBot=zeros(nt,1);
Eta=zeros(nt,1);
ThetaBot=zeros(nt,1);
Theta=zeros(nt,nz);
for t=1:nt, mydisp(t)
    fnm=['PhiBot/PhiBot_3662x392.' datestr(dt(t),30)];
    PhiBot(t)=read_ijk(fnm,nx/2,ny/2,1,nx,ny,1);
    fnm=['Eta/Eta_3662x392.' datestr(dt(t),30)];
    Eta(t)=read_ijk(fnm,nx/2,ny/2,1,nx,ny,1);
    fnm=['ThetaBot/ThetaBot_3662x392.' datestr(dt(t),30)];
    ThetaBot(t)=read_ijk(fnm,nx/2,ny/2,1,nx,ny,1);
    fnm=['Theta/Theta_3662x392x90.' datestr(dt(t),30)];
    Theta(t,:)=read_ijk(fnm,nx/2,ny/2,1:nz,nx,ny,nz);
end

region_name='Kuroshio';
pin=['~dmenemen/llc_4320/regions/' region_name '/'];
fnm=[pin 'grid/hFacC_3662x392x90'];
hFacC=read_ijk(fnm,nx/2,ny/2,1:nz,nx,ny,nz);
load ../WPac/grid/thk90
hFacC=hFacC.*thk90;
ThetaInt=0*Eta;
for t=1:nt
    ThetaInt(t)=sum(Theta(t,:).*hFacC')/sum(hFacC);
end

clear a* dp* f* p* r* nx ny* nz t h*
save figs/TimSer

PhiBot(find(PhiBot==0))=nan;

clf
subplot(211)
plot((1:nt)/24,(PhiBot-nanmean(PhiBot))/9.81*100,(1:nt)/24,(Eta-mean(Eta))*100)
axis([0 nt/24 -55 55])
legend('Bottom pressure','Sea surface height',0)
ylabel('cm'), xlabel('days from November 1, 2011')
subplot(212)
plot((1:nt)/24,(ThetaBot-mean(ThetaBot))*100,(1:nt)/24,(ThetaInt-mean(ThetaInt)))
axis([0 nt/24 -.37 .37])
legend('Bottom temperature * 100','Depth-averaged temperature',0)
ylabel('^oC'), xlabel('days from November 1, 2011')
print -dpsc TimSer

%%%%%%%%%%%%
% compare to DART data locations
% www.ngdc.noaa.gov/nndc/struts/results?&t=102597&s=1&d=1
% 46401 46.6369N -170.7976E 5540m
% 21418 38.7066N 148.6945E 5665m
% 52402 11.5687N 154.5839E 5799m
clear all, closea all
cd ~dmenemen/llc/llc_4320/regions/Kuroshio/figs
nx=4320;
kx=1:90;
prec='real*4';
gdir='~dmenemen/llc_4320/grid/';
pin='~dmenemen/llc_4320/MITgcm/run/';
fn1=[gdir 'Depth.data'];
fn2=[gdir 'hFacC.data'];
ts=10368:144:485568;
format bank

% extract indices for 46401
lat=46.6369;
lon=-170.7976;
[tmp fc ix jx]=quikread_llc(fn1,nx,1,prec,gdir,lat-.01,lat+.01,lon-.01,lon+.01);
disp([tmp fc ix jx])
% 5501.41 4 2266 10667
hFacC=read_llc_fkij(fn2,nx,fc,kx,ix,jx);
iz=max(find(hFacC));
disp([iz hFacC(iz)])
% 87 .51
ThetaBot46401=zeros(length(ts),1);
for t=1:length(ts)
    mydisp(t)
    fin=[pin 'Theta.' myint2str(ts(t),10) '.data'];
    ThetaBot46401(t)=read_llc_fkij(fin,nx,fc,iz,ix,jx);
end
clear I* a* fn* fin g* kx nx p* t tmp
save ThetaBot46401

% extract indices for 21418
lat=38.7066;
lon=148.6945;
[tmp fc ix jx]=quikread_llc(fn1,nx,1,prec,gdir,lat-.01,lat+.01,lon-.01,lon+.01);
disp([tmp fc ix jx])
% 5681.69 4 322 10109
hFacC=read_llc_fkij(fn2,nx,fc,kx,ix,jx);
iz=max(find(hFacC));
disp([iz hFacC(iz)])
% 87 1
ThetaBot21418=zeros(length(ts),1);
for t=1:length(ts)
    mydisp(t)
    fin=[pin 'Theta.' myint2str(ts(t),10) '.data'];
    ThetaBot21418(t)=read_llc_fkij(fin,nx,fc,iz,ix,jx);
end
clear I* a* fn* fin g* kx nx p* t tmp
save ThetaBot21418

% extract indices for 52402
lat=11.5687;
lon=154.5839;
[tmp fc ix jx]=quikread_llc(fn1,nx,1,prec,gdir,lat-.01,lat+.01,lon-.01,lon+.01);
disp([tmp fc ix jx])
% 5833.49 4 605 8525
hFacC=read_llc_fkij(fn2,nx,fc,kx,ix,jx);
iz=max(find(hFacC));
disp([iz hFacC(iz)])
% 88 .38
ThetaBot52402=zeros(length(ts),1);
for t=1:length(ts)
    mydisp(t)
    fin=[pin 'Theta.' myint2str(ts(t),10) '.data'];
    ThetaBot52402(t)=read_llc_fkij(fin,nx,fc,iz,ix,jx);
end
clear I* a* fn* fin g* kx nx p* t tmp
save ThetaBot52402

%%%%%%%%%%%%
clear all, close all
orient tall, wysiwyg
cd ~dmenemen/llc/llc_4320/regions/Kuroshio/figs
load ThetaBot46401
ThetaBot46401=insitutemp(33,ThetaBot46401,0,5540);
load ThetaBot21418
ThetaBot21418=insitutemp(33,ThetaBot21418,0,5665);
load ThetaBot52402
ThetaBot52402=insitutemp(33,ThetaBot52402,0,5799);
dy=(1:length(ts))/24;
subplot(311)
plot(dy,ThetaBot46401,'k')
title('DART46401 46.6369N -170.7976E 5540m (187m thick bottom layer)')
xlabel('days from September 10, 2011')
ylabel('temperature (^oC)')
subplot(312)
plot(dy,ThetaBot21418,'k')
title('DART21418 38.7066N 148.6945E 5665m (366m thick bottom layer)')
xlabel('days from September 10, 2011')
ylabel('temperature (^oC)')
subplot(313)
plot(dy,ThetaBot52402,'k')
title('DART52402 11.5687N 154.5839E 5799m (152m thick bottom layer)')
xlabel('days from September 10, 2011')
ylabel('temperature (^oC)')
print -dps DART


%%%%%%%%%%%%%

clear all, close all
cd ~dmenemen/llc_4320/regions/Kuroshio
nx=3662; ny=392; nz=90;
Tsum=0; Tssq=0; N=0;
for ts=10368:144:485568
    dy=ts2dte(ts,25,2011,9,10,30);
    fnm=['ThetaBot/ThetaBot_3662x392.' dy];
    if exist(fnm)
        disp(dy)x1g
        fld=readbin(fnm,[nx ny]);
        Tsum=Tsum+fld;
        Tssq=Tssq+fld.^2;
        N=N+1;
    end
end
clear dy f* nz ts
save figs/ThetaBotVar



%%%%%%%%%%%%%%%

clear all, close all
cd ~dmenemen/llc/llc_4320/regions/Kuroshio/figs
load NovDec
load ThetaBotVar
orient landscape
wysiwyg

clf
subplot(411)
mypcolor(Tsum'/N);
set(gca,'xticklabel',[],'yticklabel',[])
caxis([0 5])
thincolorbar
title('September-December 2011 bottom temperature (^oC)')

subplot(412)
mypcolor((ThetaBotDec'-ThetaBotNov'));
set(gca,'xticklabel',[],'yticklabel',[])
caxis([-1 1]/80)
thincolorbar
title('November to December 2011 bottom temperature change (^oC)')

subplot(413)
tmp=Tssq/N-(Tsum/N).^2;
tmp(find(tmp<0))=0;
mypcolor(sqrt(tmp'));
set(gca,'xticklabel',[],'yticklabel',[])
caxis([0 1]/40)
thincolorbar
title('Bottom temperature standard deviation (^oC)')

region_name='Kuroshio';
pin=['~dmenemen/llc_4320/regions/' region_name '/'];
fnm=[pin 'grid/Depth_3662x392'];
Depth=readbin(fnm,[3662 392]);
subplot(414)
mypcolor(-Depth'/1e3);
caxis([-6 -.5])
thincolorbar
set(gca,'xticklabel',[],'yticklabel',[])
title('Depth (km)')

print -dpng ThetaBot