Global Index (short | long) | Local contents | Local Index (short | long)
ctlim = [110 300 -30 30]; ctlim = [110 300 -72.5 67.5];
| This script calls | |
|---|---|
clear
for index = 1:3
cd /home/disk/hayes2/dvimont/csiro/data
if index == 1; tim = [101:550];
elseif index == 2; tim = [551:1000];
else; tim = [101:1000];
end
filin = ['temp_A_L1-10.nc'];
nc = netcdf(filin, 'nowrite');
depth = nc{'depth'}(:);
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
ctlim = [180 270 -6 6];
[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;
temp = squeeze(temp);
% Average data if necessary:
if 1
temp = shiftdim(temp, 1);
lat = mean([lat(1:2:nlat)'; lat(2:2:nlat)'])';
lon = mean([lon(1:2:nlon)'; lon(2:2:nlon)'])';
for i = 1:(nlat/2);
for j = 1:(nlon/2);
yk = 2*(i-1)+[1:2]; xk = 2*(j-1)+[1:2];
temp2(i,j,:) = mean2(mean2(temp(yk,xk,:)));
end
end
temp = shiftdim(temp2, 2); clear temp2;
[ntim, nlat, nlon] = size(temp);
end
temp = reshape(temp, ntim, nlat*nlon);
temp = detrend(temp);
%temp = temp ./ (ones(ntim, 1) * std(temp));
temp = cosweight(temp, lat);
clim = mean(temp);
kp = find(~isnan(clim)); nkp = length(kp);
temp = temp(:, kp);
[lam, per, lds, pcs] = eof_dan(temp);
figure(2); figure_orient;
get_global; default_global; FRAME = ctlim;
cint = 0.1; clev = -1:cint:1;
num = 1;
tem = NaN*ones(nlat*nlon, 1);
tem(kp) = lds(:,num);
tem = reshape(tem, nlat, nlon);
sd(1)
mcont(tem, clev, 'giso')
title(['EOF' num2str(num) ' of level ' num2str(num) ...
' TEMP, ' num2str(round(per(num))) '% Variance explained, ' ...
'Years ', num2str(tim(1)), ' to ', num2str(tim(ntim))]);
xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']);
sd(2)
plot(tim, pcs(:,num));
title(['PC' num2str(num)])
xlabel('Year');
ylabel('STD');
axis([min(tim)-1 max(tim)+1 -3 3]);
grid on;
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2
figure(2);% figure_orient;
for m = 1:4
subplot(4,1,m)
tind = (ntim/4)*(m-1)+[1:(ntim/4)];
plot(tind, real(pcs(tind, num)), '-k')
axis([min(tind-1) max(tind+1) -3 3]);
set(gca, 'YTick', [-3:3], 'XTick', [0:25:ntim]);
grid
end
subplot(4,1,1)
title(['LEVEL 1 TEMP: PC', num2str(num)]);
subplot(4,1,4)
xlabel('Years')
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2
index
if index == 1; pc101 = pcs;
elseif index == 2; pc551 = pcs;
end
end
% Calculate power spectrum for pcs
ctstar = detrend(squeeze(mean2(mean2(shiftdim(temp, 1)))));
ctstar = ctstar./std(ctstar);
num = 1
for i = 1:3;
if i == 1;
% ctann = pc101(:,num);
ctann = ctstar(1:450);
tit = ['Years 101 to 550'];
elseif i == 2;
% ctann = pc551(:,num);
ctann = ctstar(451:900);
tit = ['Years 551 to 1000'];
elseif i == 3;
% ctann = pcs(:,num);
ctann = ctstar(1:900);
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(1)
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 .05 15])
grid
title(['Power Spectrum for PC' num2str(num) ', 30s to 30n, ' tit]);
% title(['Power Spectrum for CT, ' tit]);
xlabel(['Frequency: yr^-^1; Degrees of Freedom: ' ...
num2str(round(dofx)) ';' ...
' Bandwidth: 7.8 x 10^-^3 yr^-^1'])
l=legend('PC1 Spectrum', 'Red Noise', '95% Confidence', '99% Confidence');
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots2