Global Index (short | long) | Local contents | Local Index (short | long)
Area average heat so process is quicker
| This script calls | |
|---|---|
clear
cd /home/disk/hayes2/dvimont/csiro/data
tim = [551:1000];
tim = [101:550];
cd /home/disk/tao/dvimont/matlab/CSIRO
ctlim = [115 285 -60 67.5];
tim = 101:550; lev = 1:7;
[temp, lat, lon, depth, middepth] = getheat(lev, tim, ctlim);
get_global; default_global; FRAME = ctlim;
[ntim, nlat, nlon] = size(temp);
heat = 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];
heat2(i,j,:) = mean2(mean2(heat(yk,xk,:)));
end
end
heat = shiftdim(heat2, 2); clear heat2;
[ntim, nlat, nlon] = size(heat);
heat = reshape(heat, ntim, nlat*nlon);
heat = detrend(heat);
heat = cosweight(heat, lat);
clim = mean(heat);
kp = find(~isnan(clim));
heat = heat(:, kp);
% Filter Heat
yr = 4.5;
[b,a]=butter(6,(2/yr));
heatlp = filtfilt(b, a, heat);
heathp = heat - heatlp;
[lam, lds, pcs, per] = complex_eof(heatlp);
lamlp = lam;
ldslp = lds;
perlp = per;
pcslp = pcs;
cd /home/disk/tao/dvimont/matlab/CSIRO/Data
% save ceof_heat1-7_yr551-1000.mat lamlp ldslp perlp pcslp kp lat lon tim ctlim
% save comp_eof_hp_heat1-7.mat lamhp ldshp perhp pcshp kp lat lon tim ctlim
% Do plots:
clear
cd /home/disk/tao/dvimont/matlab/CSIRO/Data
load comp_eof_hp_heat1-7.mat
load comp_eof_lp_heat1-7.mat
load ceof_heat1-7_yr551-1000.mat
load ceof_heat1-7_yr101-550.mat
default_global; FRAME = ctlim;
pc10 = pcs;
ld10 = lds;
ntim = size(pc10, 1);
npt = size(ld10, 1);
nlat = length(lat); nlon = length(lon);
weight1 = std(real(pc10));
weight2 = std(imag(pc10));
pc10 = real(pc10).*(ones(ntim, 1)*(1./weight1)) + ...
sqrt(-1)*imag(pc10).*(ones(ntim, 1)*(1./weight2));
ld10 = real(ld10).*(ones(npt, 1)*weight1) + ...
sqrt(-1)*imag(ld10).*(ones(npt, 1)*weight2);
figure(1); figure_landscape;
cint = 10; clims = [-100:cint:100];
num = 1;
tem = NaN * ones(1, nlat*nlon);
tem(kp) = real(ld10(:,num))';
tem = reshape(tem, nlat, nlon);
sp2(1);
% clma
gcont(tem, clims); dc;%, 'giso');
% dcm2
title(['HEAT to 275m: REAL CEOF' num2str(num)])
xlabel(['Contour Interval: ' num2str(cint) ' C (std)^-^1'])
tem = NaN * ones(1, nlat*nlon);
tem(kp) = imag(ld10(:,num))';
tem = reshape(tem, nlat, nlon);
sp2(2);
% clma
gcont(tem, clims); dc;%, 'giso');
% dcm2
title(['HEAT to 275m: IMAG CEOF' num2str(num)])
xlabel(['Contour Interval: ' num2str(4.2*cint) 'x10^6 J m^-^2 std^-^1'])
figure(2);% figure_orient;
for m = 1:4
subplot(4,1,m)
tind = [floor((ntim/4)*(m-1)+[1:(ntim/4)]) 112*m+1];
plot(tind, real(pc10(tind, num)), '-k', ...
tind, imag(pc10(tind, num)), '--k');
axis([min(tind-1) max(tind+1) -3 3]);
set(gca, 'YTick', [-3:3], 'XTick', [0:10:ntim]);
grid
end
subplot(4,1,1)
title(['AVERAGE TEMP to 620m; CPC', num2str(num), ...
': Dashed = IMAG, Solid = REAL']);
subplot(4,1,4)
xlabel('Years')
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots
% Look at various phases of this pattern
figure(1); clf;
lag = 0; lg = lag*pi/180;
num = 1;
nfrm = 6;
j = sqrt(-1);
for i = 1:nfrm;
if i == 1; figure(1); figure_landscape;
elseif i == 3; figure(2); figure_landscape;
elseif i == 5; figure(3); figure_landscape;
end
wgt = conj(exp(j * ((i-1) * pi/nfrm + lg) ) * ones(npt, 1));
pat = real(wgt .* ld10(:,num));
tem = NaN * ones(nlat*nlon, 1);
tem(kp, 1) = pat;
figure(1); clf;
lag = 0; lg = lag*pi/180;
num = 1;
nfrm = 6;
j = sqrt(-1);
for i = 1:nfrm;
if i == 1; figure(1); figure_landscape;
elseif i == 3; figure(2); figure_landscape;
elseif i == 5; figure(3); figure_landscape;
end
wgt = conj(exp(j * ((i-1) * pi/nfrm + lg) ) * ones(npt, 1));
pat = real(wgt .* ld10(:,num));
tem = NaN * ones(nlat*nlon, 1);
tem(kp, 1) = pat;
tem = reshape(tem, nlat, nlon);
sp2(2-rem(i,2));
if ismap(gca); clma; end
gcont(tem, [-200:10:200]);dc;%, 'giso');
title(['4.5YR LP CEOF' num2str(num)...
', Phase = ' num2str(((i-1)*180)/nfrm+lag) ', YEARS 101-550'])
xlabel(['Contour Interval: 10 K std^-^1'])
if ~rem(i,2)
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots3
eval(['print -dps2 LP_CEOF' num2str(num)...
'_Phase_' num2str((i-1)*180/nfrm) 'yr101-550.ps'])
end
end
for i = 3:-1:1; figure(i); end;
tem = reshape(tem, nlat, nlon);
sp2(2-rem(i,2));
if ismap(gca); clma; end
gcont(tem, [-200:10:200]);dc;%, 'giso');
title(['4.5YR LP CEOF' num2str(num)...
', Phase = ' num2str(((i-1)*180)/nfrm+lag) ', YEARS 101-550'])
xlabel(['Contour Interval: 10 K std^-^1'])
if ~rem(i,2)
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots3
eval(['print -dps2 LP_CEOF' num2str(num)...
'_Phase_' num2str((i-1)*180/nfrm) 'yr101-550.ps'])
end
end
for i = 3:-1:1; figure(i); end;
% Look at lagged correlations between real and imag cpc's
figure(4); figure_orient;
for i = 1:6
if i == 1;
v1 = real(pc10(:,1)); v2 = imag(pc10(:,1)); lab = ['Re(CPC1)'; 'Im(CPC1)'];
elseif i == 6;
v1 = real(pc10(:,1)); v2 = real(pc10(:,2)); lab = ['Re(CPC1)'; 'Re(CPC2)'];
elseif i == 3;
v1 = real(pc10(:,1)); v2 = imag(pc10(:,2)); lab = ['Re(CPC1)'; 'Im(CPC2)'];
elseif i == 4;
v1 = imag(pc10(:,1)); v2 = real(pc10(:,2)); lab = ['Im(CPC1)'; 'Re(CPC2)'];
elseif i == 5;
v1 = imag(pc10(:,1)); v2 = imag(pc10(:,2)); lab = ['Im(CPC1)'; 'Im(CPC2)'];
elseif i == 2;
v1 = real(pc10(:,2)); v2 = imag(pc10(:,2)); lab = ['Re(CPC2)'; 'Im(CPC2)'];
end
a = zeros(21,1);
for lag = -10:10;
a(lag+11) = corr(v1, v2, lag);
end
subplot(3,2,i);
bar(-10:10, a);
title(['< ' lab(1,:) ', ' lab(2,:) ' >'])
xlabel([lab(1,:) ' leads ' lab(2,:) ' | ' lab(2,:) ' leads ' lab(1,:)])
axis([-11 11 -1 1])
set(gca, 'XTick', -10:5:10, 'YTick', -1:.25:1)
grid on
end
% Look at power spectrum of CPCs:
for i = 1:2;
if i == 1; ctann = real(pc10(1:2:ntim,num)); tit = ['Real CPC' num2str(num)];
else; ctann = imag(pc10(1:2:ntim,num)); tit = ['Imag CPC' num2str(num)];
end
nx = length(ctann);
nfft = 64;
noverlap = nfft/2;
[p, f] = spectrum(ctann, nfft, noverlap);
n2 = nfft/2;
f2 = 0.25 * (0:n2)/n2;
%f2 = 2 * (0:n2)/n2;
rho = (corr(ctann, ctann, 1));% + sqrt(corr(ctann, ctann, 2)))/2;
% (corr(ctann, ctann, 3)^(1/3))) / 3;
%rho = corr(ctann, ctann, 1);
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.03;
rner5 = rn * 1.66;
dofx = nx / n2;
figure(2)
sp(i)
semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ...
f2, rner5, 'b--', f2, rner1, 'b-.')
axis([0 .25 .1 10])
grid
title(['Power Spectrum for ' tit ': SSTA(180:95w 6s:6n)'])
xlabel(['Frequency: yr^-^1; Degrees of Freedom: 15.6;'...
' Bandwidth: 7.8 x 10^-^3 yr^-^1'])
legend('CT Spectrum', 'Red Noise', '95% Confidence', '99% Confidence');
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots