%
% phase_b.m           % filter 2-d field: alongboundary vorticity vs time at 6 model 
% levels. Source file is dh3a15b5.vor referenced total vorticity)- output 
% from  benchmark run with daily forcing using <tpersf_13.dat>, 2151 days long.
% filter is "average", or "gaussian". Output file is <dh3v1015.dat>.
% ADD PATH     ADD PATH     ADD PATH       ADD PATH
       addpath /usr/people/valery/TPD_1.0/qgsbl/exam
% ***************************************************
       te=2151
       load dh3a15b5.vor
       a=dh3a15b5;
       clear dh3a15b5;
       a=a/2.3162; 
       tk=te*6;
       vort=ones(tk,108);

       for j=[1:6]
       h=ones(te,162);
       tn=j;

%       figure(j)
       p=0;
       for i=[tn:6:tk]
       p=p+1;
       a1=a(i,:);
       a1=reshape(a1,31,25);
       a1=a1';

       b=a1(1,1:31);
       h(p,24:54)=b(1:31);
       h(p,132:162)=b(1:31);
       d=a1(2:1:24,31);
       d=d';
       h(p,55:77)=d(1:1:23);

       e=a1(25,31:-1:1);
       h(p,78:108)=e(1:31);

       g=a1(24:-1:2,1);
       g=g';
       h(p,109:131)=g(1:23);
       h(p,1:23)=g(1:23);
       end
       hh=fspecial('average',[10 15]);
       h=h';
       size(h)
       hf=filter2(hh,h);
       h(132:162,:)=[];
       h(1:23,:)=[];
       hf(132:162,:)=[];
       hf(1:23,:)=[];

  dd1 = subplot(2,1,1)
       
%       contourf(300:2450,1:108,h,10);
       set(gca,'XTick',[366 731 1096 1461 1827 2192])
       set(gca,'XTickLabel',[1993 1994 1995 1996 1997 1998])
       ylabel('Boundary grid point number');
       title('Unfiltered Vorticity (sec^-^1)                       h = 100m');
       axis([300 2450 1 108]);
%      colormap copper
%       shading interp
       caxis([-50 30]);
%       colorbar

       x=300:2450;
       y1(1:te)=31;
       y2(1:te)=54;
       y3(1:te)=85;
       hold on;
%       plot(x,y1,'w-','LineWidth',1);
%       plot(x,y2,'w-','LineWidth',1);
%       plot(x,y3,'w-','LineWidth',1);
       hold off
%     
       
       subplot(2,1,2)
%       contourf(300:2450,1:108,hf,10);
       hold on;
%       plot(x,y1,'w-','LineWidth',1);
%       plot(x,y2,'w-','LineWidth',1);
%       plot(x,y3,'w-','LineWidth',1);
       hold off
       set(gca,'XTick',[366 731 1096 1461 1827 2192])
       set(gca,'XTickLabel',[1993 1994 1995 1996 1997 1998])
       title('Filtered Vorticity (sec^-^1)');
       axis([300 2450 1 108]);
       ylabel('Boundary grid point number');
       xlabel('Time (years)','fontweigh','bold');
       caxis([-50 30]);
   colormap(pink(9))
   hs=    colorbar('horiz')
   set(hs,'Position',[0.2328    0.0640    0.5750    0.0158])
   set(dd1,'Position',[0.1300    0.5811    0.7750    0.2341])
%      get back alongboundary vorticity

       hf=hf';
       p=0;
       for i=[tn:6:tk]
       p=p+1;
       vort(i,:)=hf(p,:);
       end



a1=ones(25,31);
p=0;
for i=[tn:6:tk]
p=p+1;
g1(1:1:23)=hf(p,86:1:108);
gg=g1';
a1(24:-1:2,1)=gg;

e1(1:1:31)=hf(p,55:1:85);
a1(25,31:-1:1)=e1;

d1(1:1:23)=hf(p,32:1:54);
dd=d1';
a1(2:1:24,31)=dd;

b1(1:1:31)=hf(p,1:1:31);
a1(1,1:1:31)=b1;
aa=a1';
aa=aa(:);
aa=aa';
a(i,:)=aa;

end
h2=ones(te,108);
p=0;
for i=[tn:6:tk]
p=p+1;
a2=a(i,:);
a2=reshape(a2,31,25);
a2=a2';
size(a2);
b2=a2(1,1:1:31);
h2(p,1:1:31)=b2(1:1:31);

d2=a2(2:1:24,31);
d2=d2';
h2(p,32:1:54)=d2(1:1:23);

e2=a2(25,31:-1:1);
h2(p,55:1:85)=e2(1:1:31);

g2=a2(24:-1:2,1);
g2=g2';
h2(p,86:1:108)=g2(1:1:23);

% get back
end
h2=h2';
size(h2)
%subplot(3,1,3)
%  contourf3(h2,10);
%axis([0 te 1 108]);

%caxis([-100 100]);
%colorbar

y1(1:te)=31;
y2(1:te)=54;
y3(1:te)=85;
%hold on;
%plot(x,y1,'w-','LineWidth',3);
%plot(x,y2,'w-','LineWidth',3);
%plot(x,y3,'w-','LineWidth',3);
%hold off

%ylabel('south             east              north                      west');
%xlabel('days');
%title('Phase diagram - along boundary  vorticity  vs time');

%  input('?')
end
 m='%7.2f';
     for s=[1:107]
     m=[m,' %7.2f'];
     end
     fid=fopen('dh3v1015.dat','w');
     fprintf(fid,[m,'\n'],vort');
     fclose(fid);
  

