% a matlab plotting program to read successive matrices from disk files
% this version used for basin model, among others. 

clear all
      
load what.dat
nx = what(1);
ny = what(2);
nz = what(3);
dx = what(5);
dy = what(6);
x1 = what(7);
y1 = what(9);
H = 1000.;


% rough guess at the long barotropic Rossby wave phase speed
 
beta = 1.9e-11;
k = 2*pi/(nx*dx);
c = beta/k^2
kmperday = c*8.64e4/1000.

x1 = 0.
y1 = 0.

x = zeros(1,nx);
y = zeros(1,ny);

x(:) = 1:nx;
y(:) = 1:ny;
for i=1:nx;
x(i) = x1 + (i-1)*dx;
x(i) = x(i)/1000.;
end;
for i=1:ny;
y(i) = y1 + (i-1)*dy;
y(i) = y(i)/1000.;
end;

    load basin10.dat;
    [m n] = size(basin10);
    time = zeros(m,1);   
    time(:) = basin10(:,1);
       

 load basin5.dat
 s = basin5/1.e6;

 [ntim, nx] = size(s)

 figure(1)
 clf reset

  ci = -20:2.:20;
c = contour(x, time, s, ci);
 clabel(c);
 xlabel('East distance, km')
 ylabel('time, days')
 title('Steamfuntion(y=L/2, t)')

 xq = [200 1200]
 yq = [25 25-1000/kmperday]
 hold on
 plot(xq,yq)
 %text(xq(2), yq(2), ' \beta k^{-2}')
 gtext(' \beta k^{-2}')


 figure(2)
 clf reset

 sp = zeros(1,nx); 

 vp = zeros(1,nx-1); vpx = vp;

 for j=1:5:ntim

 for k=1:nx
 sp(k) = s(j,k);
 end

 for k=1:nx-2
 vp(k) = (s(j,k+1) - s(j,k))/(H*dx);
 vpx(k) = (x(k+1) + x(k))/2.;
 end

 subplot(2,1,1)
 hold on
 plot(x, sp)
 grid
 xlabel('East distance, km')
 ylabel('Streamfunction Sv')

 subplot(2,1,2)
 hold on
 plot(vpx, vp)
 grid
 xlabel('East distance, km')
 ylabel('North velocity m s^{-1}')

 end


%  load and plot the zonal wind stress

load basin6.dat

y(:) = basin6(:,1);
x(:) = basin6(:,2)*1023.;

figure(3)
clf reset
plot(x,y)
xlabel('east wind stress')
ylabel('north distance, km')