% geoadj_2d_plot.m
% a matlab plotting program to read successive matrices from disk files
% this version used for geoadj-2d.for model. 
%
%  by Jim Price, 2003, 2010.
%

clear all
close all

path(path, 'c:\matlabextras')  %  a directory that holds mpgwrite.m (highly recommended)

%  Set some default graphics things.
set(0,'DefaultTextFontSize',18)
set(0,'DefaultAxesFontSize',18)
set(0,'DefaultAxesLineWidth',1.5)
set(0,'DefaultLineLineWidth',1.0)


% to make a decent 4 panel of eta, need to use smaller fonts......
set(0,'DefaultTextFontSize',12)
set(0,'DefaultAxesFontSize',12)

load ga_what.dat
what = ga_what;
nx = what(1);
ny = what(2);
nz = what(3);
dx = what(5);
dy = what(6);
x1 = what(7);
y1 = what(9); 
f0 = what(11);
beta = what(12);
H = what(15);

kc = 0;

x1 = 0.
y1 = 0.

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

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;

x = x - mean(x);
y = y - mean(y);
       
load ga_time.dat
[ntim, kkk] = size(ga_time);

etax = zeros(ntim, nx);
vorx = etax; 
betav = etax;

h1 = zeros(ny, nx); h2 = h1; h3 = h1;        

load ga_eta.dat;
load ga_u.dat;
load ga_v.dat; 
load ga_t.dat;

nn = 0;
nm = 0;

ii = 0;
figure(91)
clf reset

% h1 is layer 1 eta, h2, h3 are U1, V1, h5 is vorticity


for i=1:41;   %  this loop is on time

 n1 = (i-1)*ny + 1;
 n2 = n1 + ny - 1;

 h1(:,:) = ga_eta(n1:n2,:);   
 h2(:,:) = ga_u(n1:n2,:); 
 h3(:,:) = ga_v(n1:n2,:);
 h5(:,:) = ga_t(n1:n2,:);
 
 if i == 1; 
     etamax = max(max(h1));
     etamin = min(min(h1));
     if abs(etamin) >= etamax;
         etamax = etamin;
     end  
end

% try rescaling velocity to give the little ones a chance
% v = (h2.^2 + h3.^2).^0.5;
% vmax = max(max(v)) + 0.000001;

%  now subsample to put the data on a domain better for plotting

xp1 = -1500.;
dx1 = 75.;
xp2 = 500.;
yp1 = -1000.;
dy1 = 75.;
yp2 = 1000.;
xp = xp1:dx1:xp2;
yp = yp1:dy1:yp2;
[xp3, yp3] = meshgrid(xp, yp);

hp2 = interp2(x,y,h2,xp3,yp3);
hp3 = interp2(x,y,h3,xp3,yp3);

%  now for eta (to add contours to velocity)
dx1 = 20.;
dy1 = 20.;
xph = xp1:dx1:xp2;
yph = yp1:dy1:yp2;
[xp3, yp3] = meshgrid(xph, yph);

hp = interp2(x,y,h1,xp3,yp3);
 
figure(1)
clf reset
hf1 = figure(1);
set(hf1,'paperposition', [0 0 680 500])
usc = 125.*etamax;   % this should be computed from the basic paramters for each case....

%  add contours of eta
hold on
[c, hc] = contourf(xp3, yp3, hp/etamax, [-0.5:0.2:1.2]);
hq = quiver(xp, yp, usc*hp2, usc*hp3, 0.);
set(hq, 'color', 'w','linewidth', 1.5);

%  try to label the max and min of eta
clear a1 a2;
maxh = max(max(hp));
[a1] = find(hp == maxh);
xx = mean(xp3(a1)); yy = mean(yp3(a1));   
if maxh >= 5; 
ht =  text(xx, yy,'H');
set(ht, 'fontsize', 12, 'fontweight', 'bold', 'color', 'w'); 
end

clear a1 a2;
minh = min(min(hp));
[a1] = find(hp == minh);
a2 = a1(1);
if minh <= -1; 
    ht = text(xp3(a2), yp3(a2),'L');
    set(ht, 'fontsize', 12, 'fontweight', 'bold', 'color', 'w'); 
end

xlabel('east, km')
ylabel('north, km')
 axis([-1500 500 -1000 1000])
 tim = ga_time(i);
 text(-800., 650., ['time = ',num2str(tim,3), ' days '], 'color', 'w');
 nf1 = 560
 nf2 = 418
 nm = nm + 1;
 Mu(nm) = getframe(hf1, [0 0 nf1 nf2]);
 nm = nm + 1;
 Mu(nm) = getframe(hf1, [0 0 nf1 nf2]);
  nm = nm + 1;
 Mu(nm) = getframe(hf1, [0 0 nf1 nf2]);

% add the equator

 xx = [0. x(nx)];
 yd2 = y(ny)/2.;
 yy = [yd2 yd2];
 hold on
% line(xx,yy)

if i == ntim
 saveas(hf1, 'ga_long_eta.fig')     %  save the last frame
end

if i == 2
    hq = h1(2,2);
end
h1 = h1 - h2;

%   now plot the tracer ...............
notr = 0
if notr == 0
hf02 = figure(20)
clf reset
set(gcf,'paperposition', [0 0 825 500])

 % pcolor(h1)
 
dx1 = 20.;
dy1 = 20.;
xp = xp1:dx1:xp2;
yp = yp1:dy1:yp2;
[xp3, yp3] = meshgrid(xp, yp);

hp = interp2(x,y,h5,xp3,yp3);
% hp = hp/500.;

  mesh(xp3, yp3, hp);

% contour(xp3, yp3, h5)
caxis([0 1])
 
 % view(-30., 50.)
 % axis([-1500 500 -1000 1000 0 1.5 ])
  
 xlabel('east, km ');  
 ylabel('north, km ');
 zlabel('q');  
 
 % add the time
 tim = ga_time(i) - 0.007;
 tim = round(10*tim)/10;
 text(-800., 600., 1.2, ['time = ',num2str(tim,3), ' days '])      
 pause(0.5)
end    %  on plot tracer 

%   now plot the displacement...............
hf2 = figure(2)
clf reset
set(gcf,'paperposition', [0 0 680 500])

 % pcolor(h1)
 
dx1 = 20.;
dy1 = 20.;
xp = xp1:dx1:xp2;
yp = yp1:dy1:yp2;
[xp3, yp3] = meshgrid(xp, yp);

hp = interp2(x,y,h1,xp3,yp3);
 
 mesh(xp3, yp3, hp/etamax);
 caxis([-0.2 1.2])
 
 view(-30., 50.)
 axis([-1500 500 -1000 1000 -0.4 1.2 ])
  
 xlabel('east, km ');  
 ylabel('north, km ');
 zlabel(['\eta/\eta_{max}']) % , \eta_{max} =' num2str(etamax,2), ' m']);           
  
 % add the time
 tim = ga_time(i);
 text(-800., 600., 1.1, ['time = ',num2str(tim,3), ' days '])                            
 pause(0.5);
 nn = nn + 1

 if (i == 1 | i == 5 | i == 13 | i == 39)
 ii = ii + 1;   
     hf9 = figure(91)
     if ii == 1
     set(hf9,'paperposition', [0 0 8 6]) 
     end
     subplot(2,2,ii) 
     mesh(xp3, yp3, hp/etamax);
     caxis([-0.4 1.2])
     view(-30., 50.)
 axis([-1500 500 -1000 1000 -0.4 1.2 ])
  
 xlabel('east, km ');  
 ylabel('north, km ');
 zlabel('\eta/\eta_{max}');           
  
 % add the time
 tim = ga_time(i) - 0.007;
 tim = round(10*tim)/10;
 text(-800., 600., 1.2, [num2str(tim,3), ' days '])      
     
 end   % if on ii
 
 
 
 nyc = round(ny/2);
 
     kc = kc + 1;
     etax(kc,:) = h1(nyc,:);
%      vorx(kc,:) = H*h5(nyc,:) - f0;
     vorx(kc,:) = h5(nyc,:);
     betav(kc,:) = h3(nyc,:);
     
 figure(25)
 clf reset
 subplot(2,1,1)
 plot(x, etax(kc,:))

 subplot(2,1,2)
 plot(x, vorx(kc,:))


 figure(2)
 nf1 = 570
 nf2 = 418
 Mh(nn) = getframe(hf2, [0 0 nf1 nf2]);
 nn = nn + 1
 Mh(nn) = getframe(hf2, [0 0 nf1 nf2]);
 nn = nn + 1
 Mh(nn) = getframe(hf2, [0 0 nf1 nf2]);

end
 
%  evaluate the terms in linear q-conservation

nc = 30
i = nc

dt7 = 8.64e4*(looktim(i+1) - looktim(i-1));
g = 9.8;
dr = 2.0;
r0 = 1030.;
eta0 = 20.;
C = sqrt(g*dr*H/r0);
bu = beta*(C*eta0/H);

dvor = (vorx(nc+1,:) - vorx(nc-1,:))/dt7;
betav5 = beta*betav(nc,:);
deta = f0*(etax(nc+1,:) - etax(nc-1,:))/(dt7*H);



  
imov = input('  Shall we save mpeg movie files?  (1 = yes)')
if imov == 1
    colorm = colormap;
if  exist('geoadj-2d-h.m') == 1; 
    delete 'geoadj-2d-h.mpg';
    delete 'geoadj-2d-u.mpg'
end
    mpgwrite(Mh, colorm, 'geoadj-2d-h.mpg')
    mpgwrite(Mu, colorm, 'geoadj-2d-u.mpg')
end