% Correlate November values of the QBO and MSU channel 4 temperature
%
% 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 ) * 0.1;

write_ts( ts*10, yr1, yr2 );

% Apparently the MSU data is only through August 1997.
[ junk, i4, r4, nam, xgrid, ygrid ] = nc_read( 'msu4.mon.anom.7997.nc' );
range( junk );
plot_image10c( xgrid, ygrid, junk ); drawcoasts

ts = time_shift( ts, yr1, yr2, 1979, 1996 );
yr1 = 1979
yr2 = 1996
nyr = yr2 - yr1 + 1
write_ts( ts*10, 1979, 1996 )

fdat = zeros(nyr,length(junk));
whos fdat
for icnt = 1: nyr
fdat(icnt,:) = nc_read( 'msu4.mon.anom.7997.nc', (icnt-1)*12+11 );
end
length(find( isnan(mean(fdat)) ))
fdat = fdat - ones(nyr,1) * mean(fdat);

tsmean = mean(ts(11:12:nyr*12))
tsstd  = std( ts(11:12:nyr*12))
ts = ( ts - tsmean ) / tsstd;
% ts = ts(11:12:nyr*12);


reg = ts(11:12:nyr*12)' * fdat / ( nyr - 1 );
stdmap = std(fdat);
range( stdmap );
cor = reg ./ stdmap;
range( cor );

plot_image10c( xgrid, ygrid, cor, 1 ); drawcoasts

i4(1:5) = [ 1 1 1 0 0 ]
r4(1:2) = [ 0 0.001 ];
nc_write( 'qbomsu4correg.nc', 1, cor, '1979-96 Novembers K', i4, r4 );
nc_write( 'qbomsu4correg.nc', 2, reg, '1979-96 Novembers K', i4, r4 );

colorvec = ones(size(ts));
colorvec(11:12:nyr*12) = colorvec(11:12:nyr*12) - ones(nyr,1);

clf
fillbar3( ts, 1979, 1996, 12, 1979, colorvec );
set( gca, 'YLim', [ -2 2 ] );
set( gca, 'DataAspectRatio', [ 100 5 1 ] )
set( gca, 'XTickLabel', [ '79'; '80'; '  '; '82'; '  '; '84'; '  '; ...
  '86'; '  '; '88'; '  '; '90'; '  '; '92'; '  '; '94'; '  '; '96'; ...
  '  '; '98'; '  '; '00' ] );
font( 'Times', 12 );
title2( 'QBO u-wind anoms., 1979-96, stand. by Novembers (black)' )