% Read in Naujokat's QBO index, make a plot of it, 
% and write out some files of this data.
%
% The first step is to take the table received from
% Barbara Naujokat and remove all of the character
% data so that it is a simple table.  Matlab wants
% all of the rows to be the same number of columns.
% Call this data edited data file "naujokat.dat".
%
% Todd Mitchell, February 1999

load naujokat.dat

yr1 = naujokat(1,1)
yr2 = naujokat(size(naujokat,1),1)
nyr = yr2 - yr1 + 1
nt  = ( yr2 - yr1 + 1 ) * 12

ts = reshape( naujokat(:,2:13)', nt, 1 );

temp = ts * 0.1;
colorvec = ones(size(temp)) * 2;
x = find( temp>0 );
colorvec(x) = colorvec(x) + ones(size(x));

clf
for iyr = 1950: yr2+1
xval = (iyr-1950)*12+1;
plot( [ xval xval ], [ -40 40 ], 'k-', 'LineWidth', 1, ...
  'Color', ones(3,1)*0.9 )
hold on
end
fillbar3( temp, yr1, yr2, 12, 1950, colorvec );
set( gca, 'YLim', [ -40 40 ] );
ytick( [ -40: 10: 40 ] );
set( gca, 'PlotBoxAspectRatio', [ 3 1 1 ] );
set( gca, 'XTicklabel', [ '1950'; '    '; '    '; '    '; '    '; ...
  ' 55 '; '    '; '    '; '    '; '    '; ...
  ' 60 '; '    '; '    '; '    '; '    '; ...
  ' 65 '; '    '; '    '; '    '; '    '; ...
  ' 70 '; '    '; '    '; '    '; '    '; ...
  ' 75 '; '    '; '    '; '    '; '    '; ...
  ' 80 '; '    '; '    '; '    '; '    '; ...
  ' 85 '; '    '; '    '; '    '; '    '; ...
  ' 90 '; '    '; '    '; '    '; '    '; ...
  ' 95 '; '    '; '    '; '    '; '    '; '2000'; '    '; '    ' ] );
font( 'Times', 12 );
title2( 'QBO 30hPa zonal wind (m/s), 1953 - Sep 2001' )
text( 1, -70,  ...
  [ 'Canton Island (3S,172W) Jan 1953 - Aug 1967 ' setstr(10) ...
    'Gan/Maledives (1S, 73E) Sep 1967 - Dec 1975' setstr(10) ...
    'Singapore     (1N,104E) Jan 1976 - Sep 2001' ], ...
  'FontName', 'Times', 'FontSize', 10 );
a = get( gca, 'XLim' );
text( a(2)-12, -65,  'Freie Universitaet Berlin', ...
  'FontName', 'Times', 'FontSize', 10, 'FontAngle', 'italic', ...
  'HorizontalAlignment',  'right' );
print -dpsc2 qbo.ps
print -djpeg qbo.jpg

[ years, months ] = yearsmonths( yr1, yr2 );
temp = ts * 0.1;
a = [ years months temp ];
save qbo19532001.asci a -ascii
save qbo19532001.dat temp -ascii