Global Index (short | long) | Local contents | Local Index (short | long)
ctlim = [-0.1 17.5 47.5 80]; ctlim = [300 360 47.5 80]; ctlim = [110 300 -30 30];
| This script calls | |
|---|---|
clear
cd /home/disk/hayes2/dvimont/csiro/data
tim = [101:550];
filin = ['temp_A_L1-10.nc'];
nc = netcdf(filin, 'nowrite');
depth = nc{'depth'}(:);
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
ctlim = [110 300 -67.5 60];
% ctlim = [110 300 12.5 62.5];
% ctlim = [110 300 -22.5 22.5];
% ctlim = [110 300 22.5 62.5];
[xk, yk] = keep_var(ctlim, lon, lat);
temp = nc{'temp'}(tim,:,yk,xk);
mv = nc{'temp'}.missing_value(:);
nc = close(nc);
lev2 = 1:3:10;
temp = temp(:,lev2,:,:);
temp(find(temp == mv)) = NaN * ones(size(find(temp == mv)));
lat = lat(yk); lon = lon(xk);
[ntim, nlev, nlat, nlon] = size(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 l = 1:nlev;
for i = 1:(nlat/2);
% for i = 1:nlat;
for j = 1:(nlon/2);
yk = 2*(i-1)+[1:2]; xk = 2*(j-1)+[1:2];
% yk = i; xk = 2*(j-1)+[1:2];
temp2(l,i,j,:) = mean2(mean2(squeeze(temp(l,yk,xk,:))));
% temp2(l,i,j,:) = mean2(squeeze(temp(l,yk,xk,:)));
end
end
end
end
temp = shiftdim(temp2, 3); clear temp2;
[ntim, nlev, nlat, nlon] = size(temp);
% Do this on the correlation matrix:
temp = reshape(temp, ntim, nlev*nlat*nlon);
for i = 1:nlev*nlat*nlon;
temp(:,i) = temp(:,i)./std(temp(:,i));
end
temp = reshape(temp, ntim, nlev, nlat, nlon);
% Weight by depth
if 0;
middepth = [0; (depth([lev2])+depth(lev2+1)) / 2];
wgt = diff(middepth)/100;
wgt = wgt./mean(wgt);
temp = shiftdim(temp, 1);
temp = reshape(temp, nlev, nlat*nlon*ntim);
temp = temp .* (wgt * ones(1, nlat*nlon*ntim));
temp = reshape(temp, nlev, nlat, nlon, ntim);
temp = shiftdim(temp, 3);
end
for i = 1:nlev;
temp(:,i,:,:) = cosweight(squeeze(temp(:,i,:,:)), lat);
end;
temp = reshape(temp, ntim, nlev*nlat*nlon);
temp = detrend(temp);
clim = mean(temp);
kp = find(~isnan(clim)); npt = length(kp);
temp = temp(:, kp);
% Lowpass filter
if 1
yr = 8;
[b,a]=butter(6,(2/yr));
temp = temp - filtfilt(b, a, temp);
end
if 0;
c = temp * temp' / npt;
[lam, pcs, per] = eof(c);
pc10 = pcs(:,1:10);
pc10 = (pc10 - ones(ntim, 1)*mean(pc10)) ./ (ones(ntim, 1)*std(pc10));
ld10 = pc10' * temp ./ ntim;
cd /home/disk/hayes2/dvimont/csiro/matlab_data
%save temp_3deof_pac.mat pc10 ld10 kp lat lon lev2 lam per
clear c pcs
end
cd /home/disk/hayes2/dvimont/csiro/matlab_data
[lam, lds, pcs, per] = complex_eof(temp);
%save 6yrLP_PAC_NATL_CEOF.ps lam lds pcs per kp lat lon lev2 depth kp kp2 kp3 lat lon lat2 lon2 lat3 lon3
%save LPtemp_3dceof_22.5to62.5pac.mat lam lds pcs per kp lat lon lev2 depth
%save HPtemp_3dceof_12.5to62.5pac.mat lam lds pcs per kp lat lon lev2 depth
%save 8yrHPtemp_3dceof_-67.5to60pac.mat lam lds pcs per kp lat lon lev2 depth
load tem_latlon.mat
%
%
%
clear
cd /home/disk/hayes2/dvimont/csiro/matlab_data
load temp_3dceof_pac.mat
get_global; default_global; FRAME = [110 300 -75 65];
ntim = size(pcs, 1);
npt = length(kp);
nlev = length(lev2);
nlat = length(lat);
nlon = length(lon);
% Weight the CPCs and CEOFs accordingly
weight1 = std(real(pcs));
weight2 = std(imag(pcs));
pcs = real(pcs).*(ones(ntim, 1)*(1./weight1)) + ...
j*imag(pcs).*(ones(ntim, 1)*(1./weight2));
lds = real(lds).*(ones(npt, 1)*weight1) + ...
j*imag(lds).*(ones(npt, 1)*weight2);
% Reshape temperature field
temp = NaN*ones(10, nlev*nlat*nlon);
temp(:,kp) = lds';
temp = reshape(temp, 10, nlev, nlat, nlon);
% Unweight the depth dimension
middepth = [0; (depth([lev2])+depth(lev2+1)) / 2];
wgt = diff(middepth)/100;
wgt = wgt./mean(wgt);
temp = shiftdim(temp, 1);
temp = reshape(temp, nlev, nlat*nlon*10);
temp = temp ./ (wgt * ones(1, nlat*nlon*10));
temp = reshape(temp, nlev, nlat, nlon, 10);
temp = shiftdim(temp, 3);
for lev = 1:5
figure(lev); figure_orient;
num = 1;
tem = real(squeeze(temp(num, lev, :, :)));
cint = 0.05; clev = [-.5:cint:.5];
sp(1)
if ismap(gca); clma; end;
mcont(tem, clev, 'giso')
title(['3D TEMP: Real CEOF' num2str(num) ', ' ...
num2str(depth(lev2(lev))/100) 'm DEPTH, ' ...
num2str(round(per(num))) '% Variance Explained']);
xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']);
tem = imag(squeeze(temp(num, lev, :, :)));
cint = 0.05; clev = [-.5:cint:.5];
sp(2)
if ismap(gca); clma; end;
mcont(tem, clev, 'giso')
title(['3D TEMP: Imag CEOF' num2str(num) ', ' ...
num2str(depth(lev2(lev))/100) 'm DEPTH, ' ...
num2str(round(per(num))) '% Variance Explained']);
xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1']);
end
for lev = 5:-1:1; figure(lev); end;
figure(6);% 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', ...
tind, imag(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(['3D TEMP: CPC', num2str(num), ...
': Dashed = IMAG, Solid = REAL']);
subplot(4,1,4)
xlabel('Years')
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots
% Out of curiosity, look at spectrum of CPC:
for i = 1:2;
if i == 1;
ctann = real(pcs(1:450,num));
elseif i == 2;
ctann = imag(pcs(1:450,num));
end
nx = length(ctann);
nfft = 100;
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));
rner1 = rn * 2.05;
rner5 = rn * 1.68;
dofx = nx / n2;
figure(1)
sp(i)
semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ...
f2, rner5, 'b--', f2, rner1, 'b-.')
set(gca, 'YTick', [.01 .1 1 10 100]);
axis([0 0.5 .05 15])
grid
title(['Power Spectrum for CPC' num2str(num) ', 3D TEMP']);
xlabel(['Frequency: yr^-^1; Degrees of Freedom: ' ...
num2str(round(dofx)) ';'...
' Bandwidth: 7.8 x 10^-^3 yr^-^1'])
l=legend('CPC Spectrum', 'Red Noise', '95% Confidence', '99% Confidence');
end