Documentation of compare_cts2

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

  temp = nc{'temp'}(tim, yk, xk);

Cross-Reference Information

Listing of script compare_cts2

clear
cd /home/disk/hayes2/dvimont/csiro/data
tim = [101:1000];
filin = ['temp_A_L1-10.nc'];
nc = netcdf(filin, 'nowrite');
  depth = nc{'depth'}(1:7);
  lat = nc{'latitude'}(:);
  lon = nc{'longitude'}(:);
  ctlim = [-.01 360 -90 90];
  [xk, yk] = keep_var(ctlim, lon, lat);
  temp = nc{'temp'}(tim,1,yk,xk);
  mv = nc{'temp'}.missing_value(:);

nc = close(nc);
[ntim, nlev, nlat, nlon] = size(temp);
temp(find(temp == mv)) = NaN * ones(size(find(temp == mv)));
lat = lat(yk); lon = lon(xk);
get_global; default_global; FRAME = ctlim;

ctlim = [180 270 -6 6];
[xk, yk] = keep_var(ctlim, lon, lat);

ct = squeeze(mean2(mean2(shiftdim(temp(:,:,yk,xk), 1))));
ct = detrend(ct);

temp = reshape(temp, ntim, nlat*nlon);
temp = detrend(temp);
clim = mean(temp);
kp = find(~isnan(clim));
temp = temp(:,kp);

tim = [101:550];
ctem = detrend(ct(tim)); ctem = ctem./std(ctem);

ctpat = ctem' * detrend(temp(tim, :)) ./ ntim;
ct2 = ctpat * detrend(temp(tim,:))' ./ (length(kp));
ct2 = ct2./std(ct2);
ctpat2 = ct2 * detrend(temp(tim, :)) ./ ntim;

get_global; default_global;
figure(1)
sp(1)
     tem = NaN * ones(1, nlat*nlon);
     tem(kp) = ctpat;
     tem = reshape(tem, nlat, nlon);
     gcont(tem, [-.5:.05:.5]);
     dc

sp(2)
     tem = NaN * ones(1, nlat*nlon);
     tem(kp) = ctpat2;
     tem = reshape(tem, nlat, nlon);
     gcont(tem, [-.5:.05:.5]);
     dc

for i = 1:3;
  if i == 1;
    ctann = pcs1(:,1)/std(pcs1(:,1));
    tit = 'PC1, Years 101 - 550';
  elseif i == 2;
    ctann = pcs2(:,1)/std(pcs2(:,1));
    tit = 'PC1, Years 551 - 1000';
  elseif i == 3;
    ctann = pcs3(:,1)/std(pcs3(:,1));
    tit = 'PC1, Years 101 - 1000';
  end

  nx = length(ctann);

  nfft = 128; 
  noverlap = nfft/2;

  [p, f] = spectrum(ctann, nfft, noverlap);
  n2 = nfft/2;
  f2 = 0.5 * (0:n2)/n2;
  rho = (corr(ctann, ctann, 1) + sqrt(corr(ctann, ctann, 2))) / 2 % + ...
%        (corr(ctann, ctann, 3)^(1/3))) / 3;
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn);
  rn = rn * mean(p(:,1));
  if i ~= 3;
    rner1 = rn * 2.64;
    rner5 = rn * 2.01;
  else
    rner1 = rn * 2.05;
    rner5 = rn * 1.68;
  end
  dofx = nx / n2; 
  figure(2); figure_orient;
  subplot(3,1,i)
     semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ...
	      f2, rner5, 'b--', f2, rner1, 'b-.')
%     tim = 2:(n2+1); f2 = f2(tim); p = p(tim,:); rn = rn(tim);
%     rner1 = rner1(tim); rner5 = rner5(tim);
%     plot(1*log(f2), f2'.*p(:,1), 'b-', 1*log(f2), f2.*rn, 'b-', ...
%	      1*log(f2), f2.*rner5, 'b--', 1*log(f2), f2.*rner1, 'b-.')
     axis([0 0.5 .05 10])
     grid
%     title(['Power Spectrum for PC' num2str(num) ', 30s to 30n, ' tit]);
%     title(['Power Spectrum for CT' num2str(fig) ', ' tit]);
     title(tit)
     xlabel(['Frequency:  yr^-^1;  Degrees of Freedom:  ' ...
              num2str(round(dofx)) ';' ...
             '  Bandwidth:  7.8 x 10^-^3 yr^-^1'])
%     l=legend('CT Spectrum', 'Red Noise', '95% Confidence', '99% Confidence');
end

cd /home/disk/tao/dvimont/matlab/CSIRO/Plots3


%  Check out EOF of entire globe's SST

temp = squeeze(shiftdim(temp, 1));

for i = 1:nlat/2
  indy = 2*(i-1)+[1:2];
  for j = 1:nlon/2
    indx = 2*(j-1)+[1:2];
    sst(:,i,j) = squeeze(mean(mean(temp(indy,indx,:))));
  end
end
lat = mean([lat(1:2:nlat)'; lat(2:2:nlat)']);
lon = mean([lon(1:2:nlon)'; lon(2:2:nlon)']);

tim = [451:900];
tim = [1:450];
tim = [1:900];
temp = sst(tim,:,:);

[ntim, nlat, nlon] = size(temp);
temp = cosweight(temp, lat);
temp = reshape(temp, ntim, nlat*nlon);
clim = mean(temp);
kp = find(~isnan(clim));
temp = temp(:,kp);
temp = detrend(temp);

[lam, per, lds, pcs] = eof_dan(temp);
%pcs1 = pcs;
%pcs2 = pcs;
pcs3 = pcs;

ctem = pcs1(:,1);
ct2 = pcs2(:,1);