%  plot some time series data generated by basin.f
%  by Jim Price.
% 
        
clear all

load basin40.dat         %  this is time series of u,v,h
[m n] = size(basin40)

n = n-7
timez = zeros(m, 1);
timez = basin40(:,1);

uz = zeros(m, 3); vz = uz; hz = uz; 
vea = zeros(m, 1); vwa = vea; curlt = vwa;

for j=1:(n-1)/3
%uz(:,j) = basin40(:,2+(j-1)*3);
vz(:,j) = basin40(:,3+(j-1)*3);
hz(:,j) = basin40(:,4+(j-1)*3);
end

vea(:) = -basin40(:,11);
vwa(:) = -basin40(:,16);
curlt(:) = basin40(:,17);

figure(1) 
clf reset

%subplot(3,1,1)
%plot(timez, vz, timez, vea, '+')
subplot(3,1,2)
plot(timez, vwa)
ylabel('wbc transport')
grid
subplot(3,1,1)
plot(timez, curlt)
ylabel('curl tau')
grid


load basin10.dat

[nt xx] = size(basin10);
time = ones(nt,1);
u = time; v = time; h = time; t = time;
rke = time; rkei = time; rpe = time;

time(:) = basin10(:,1);
dt = time(2) - time(1)

u(:) = basin10(:,2);
v(:) = basin10(:,3);
h(:) = basin10(:,4);
t(:) = basin10(:,5);

figure(2)
clf reset

subplot(3,1,1)
plot(time,u, time,v)
xlabel('time, days')
ylabel('currents, m/sec')

subplot(3,1,2)
plot(time,t)
ylabel('temperature')
             
rke(:) = basin10(:,6);
rkei(:) = basin10(:,7);
rpe(:) = basin10(:,8);

subplot(3,1,3)
plot(time,rke,'-', time,rkei,':', time,rpe,'-.') 
xlabel('time, days')
ylabel('kinetic energy') 




ah = time; ww = time; fw = time; bw = time;
enet = time;

% plot average layer depth

figure(3)
clf reset

ah(:) = basin10(:,11);
plot(time, ah)
xlabel('time, days')
ylabel('average layer thickness, m')
 



figure(4)
clf reset
          
enet(:) = basin10(:,9); 
ww(:) = basin10(:,10);
fw(:) = basin10(:,11);
bw(:) = basin10(:,12);

plot(time, ww, time, fw, time, bw)
legend('wind work', 'friction', 'bottom')


figure(5)
clf reset

fs = time(6) - time(5)

[px, f] = psd(u, 256, fs, 256, 128);
[py, f] = psd(v, 256, fs, 256, 128);

subplot(2,1,1)
loglog(f, px, f, py, ':')
ylabel('spectra, MKS')
title('current spectra')

%  plot the frequency of the lowest basin Ro mode  
r1 = 32.58
r1f = 1/r1
xf(1) = r1f; xf(2) = r1f; yf(1) = 1.; yf(2) = 1.e-20;
hold on
loglog(xf, yf, ':')

subplot(2,1,2)
%[hx, f] = psd(h, 256, fs, 256, 128);
[tx, f] = psd(t, 256, fs, 256, 128);
loglog(f, tx, ':')
title('temperature spectrum')
xlabel('frequency, 1/days')
ylabel('spectra, MKS')



figure(6)
clf reset

[tv, f] = csd(t, v, 256, fs, 256, 128);
[ctv, f] = cohere(t, v, 256, fs, 256, 128);
tvr = abs(real(tv));
tvmag = (abs(tv.^2).^0.5);
subplot(3,1,1)
loglog(f, tvr, f, tvmag, ':')
title(' t v cross-spectra')
subplot(3,1,2)
semilogx(f,ctv);
xlabel('frequency, 1/days')
ylabel('coherence')

subplot(3,1,3)
deltf = 1./(time(256) - time(1));
tvcum = cumsum(tv)*deltf;
semilogx(f, tvcum)
ylabel('t*v sum')





load basin40.dat
[m n] = size(basin40)

timez = ones(m, 1);
timez(:) = basin40(:,1);

uz = ones(m, n);
vz = ones(m, n);

for j=1:(n-1)/3
uz(:,j) = basin40(:,2+(j-1)*3);
vz(:,j) = basin40(:,3+(j-1)*3);
end

figure(7) 
clf reset

subplot(2,1,1)
plot(timez, uz)
subplot(2,1,2)
plot(timez, vz)


figure(8)
clf reset

[px, f] = psd(u, 1028, fs, 1028, 100);
[py, f] = psd(v, 1028, fs, 1028, 100);

loglog(f, px, f, py, '--')
axis([1.e-2 5.e-1, 10e-11, 10e-4])
ylabel('spectra, MKS')
title('current spectra')

%  plot the frequency of the lowest (basinwide) Ro modes  


beta = .198e-10

xx = 1.e6; yy = 1.e6


for jj = 1:3
for kk = 1:3

bsize =  sqrt( 1/(jj*xx)^2 + 1/(kk*yy)^2 )
per = ((2*pi)^2*bsize/beta)/8.64e4
freq = 1/per

xf(1) = freq; xf(2) = freq; yf(1) = 1.; yf(2) = 1.e-20;
hold on
loglog(xf, yf, ':')

end   %  loop on jj
end   %  loop on kk

cortv = corrcoef(t, v)
covtv = cov(t, v)