Documentation of compare_cts

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

Help text

tim = [1200:12000];
filin = ['temp_M_L1_1000_years_new.nc'];

Cross-Reference Information

Listing of script compare_cts

clear
cd /home/disk/hayes2/dvimont/csiro/data
tim = [101:1000];

filin = ['temp_A_L1-10.nc'];
for ind = 1:3;
nc = netcdf(filin, 'nowrite');

  depth = nc{'depth'}(1:7);
  lat = nc{'latitude'}(:);
  lon = nc{'longitude'}(:);

  if ind == 1;
    ctlim = [180 270 -6 6];
  elseif ind == 2;
    ctlim = [165 225 -8 8];
  elseif ind == 3;
    ctlim = [175 200 25 40];
  end

  [xk, yk] = keep_var(ctlim, lon, lat);

%  temp = nc{'temp'}(tim, yk, xk);
  temp = nc{'temp'}(tim,1,yk,xk);
  mv = nc{'temp'}.missing_value(:);

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

temp = squeeze(temp);

eval(['ct' num2str(ind) ' = squeeze(mean2(mean2(shiftdim(temp, 1))));']);
eval(['[ct' num2str(ind) ', clim' num2str(ind) ...
      '] = annave(ct' num2str(ind) ');']);
end

if 0;
for i = 1:(ntim-1)/3;
  ind = 3*(i-1) + [1:3];
  ct1a(i) = mean(ct1(ind));
  ct2a(i) = mean(ct2(ind));
  ct3a(i) = mean(ct3(ind));
end
ct1 = ct1a;
ct2 = ct2a;
ct3 = ct3a;
ntim = length(ct1);
end;

for fig = 1:3;
eval(['ctstar = ct' num2str(fig) ';']);

for i = 1:3;
  if i == 1;
    ctann = ctstar(1:ntim/2);
%    ctann = ctstar(ntim/2:3*ntim/4);
    tit = ['Years 101 to 550'];
  elseif i == 2;
    ctann = ctstar(ntim/2:ntim);
%    ctann = ctstar((3*ntim/4+1):ntim);
    tit = ['Years 551 to 1000'];
  elseif i == 3;
    ctann = ctstar(1:ntim);
%    ctann = ctstar(ntim/2:ntim);
    tit = ['Years 101 to 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(fig); 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-.')
     set(gca, 'YTick', [.01 .1 1 10 100]);
     axis([0 0.5 .005 5])
     grid
%     title(['Power Spectrum for PC' num2str(num) ', 30s to 30n, ' tit]);
     title(['Power Spectrum for CT' num2str(fig) ', ' 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

end
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2

for i = 1:3;
  eval(['ct' num2str(i) ' = detrend(ct' num2str(i) ');']);
end

ct = (ct1-mean(ct1));
ct = (ct2-mean(ct2));
ct = (ct3-mean(ct3));

figure(4); figure_orient;
for i = 1:4;
  ind = 250*(i-1)+[1:250];
  subplot(4,1,i)
  plot(ind, ct(ind));
  grid on;
  axis([min(ind-1) max(ind+1) -1 1]);
end

for i = 1:900;
  ind = (i+50)+[-50:50];
  meanct(i) = mean(ct(ind));
  stdct1(i) = std(ct(ind));
end
figure(5);% figure_landscape;
sp(1)
  plot(51:950, meanct);
sp(2)
  plot(51:950, stdct1);
  axis([0 1001 0 0.35])