Global Index (short | long) | Local contents | Local Index (short | long)
tim = 551:1000;
| This script calls | |
|---|---|
clear
cd /home/disk/tao/dvimont/matlab/CSIRO/Heat
ctlim = [110 300 -60 -20];
tim = 101:550;
[heat, lat, lon, depth, middepth] = getheat(1:7, tim, ctlim);
[ntim, nlat, nlon] = size(heat);
% Possibly average the data:
if 0;
heat2 = NaN * ones(ntim, nlat/2, nlon/2);
for i = 1:nlat/2;
latind = 2*(i-1) + [1:2];
for j = 1:nlon/2;
lonind = 2*(j-1) + [1:2];
heat2(:,i,j) = mean2(mean2(shiftdim(heat(:,latind,lonind), 1)));
end
end
lat2 = mean([lat(1:2:nlat)'; lat(2:2:nlat)'])';
lon2 = mean([lon(1:2:nlon)'; lon(2:2:nlon)'])';
end
%for biff = 1:3;
biff = 2;
if biff == 1;
heat2 = heat;
elseif biff == 2;
[b, a] = butter(6, 2/10);
[b2, a2] = butter(6, 2/60);
heat2 = filtfilt(b, a, heat);% - filtfilt(b2, a2, heat);
elseif biff == 3;
[b, a] = butter(6, 2/8);
heat2 = heat - filtfilt(b, a, heat);
end
lon2 = lon;
lat2 = lat;
[ntim, nlat, nlon] = size(heat2);
[heat2, clim] = remove_mean(heat2);
heat2 = reshape(heat2, ntim, nlat*nlon);
heat2 = detrend(heat2);
heat2 = reshape(heat2, ntim, nlat, nlon);
% Normalize
if 0;
heat2 = reshape(heat2, ntim, nlat*nlon);
heat2 = heat2 ./ (ones(ntim,1) * std(heat2));
heat2 = reshape(heat2, ntim, nlat, nlon);
end
heat2 = cosweight(heat2, lat2);
kp = find(~isnan(clim));
heat2 = heat2(:,kp);
%[lam, lds, pcs, per] = complex_eof(heat2);
[lam, lds, pcs, per] = eof_dan(heat2);
disp('done!');
cd /home/disk/hayes2/dvimont/csiro/matlab_data/Heat_Content
if biff == 1;
save RAW_L1-7_EOF.mat lam lds pcs per ctlim lat lon
elseif biff == 2;
save SP2_LP10_detrend_L1-7_yr101-550.mat lam lds pcs per ctlim lat lon
elseif biff == 3;
save HP8_L1-7_EOF.mat lam lds pcs per ctlim lat lon
end
end
% Look at CEOFs
clear
cd /home/disk/tao/dvimont/matlab/CSIRO/Heat
ctlim = [110 300 -60 60];
[heat, lat, lon, depth, middepth] = getheat(1:7, 101:550, ctlim);
[ntim, nlat, nlon] = size(heat);
heat = reshape(heat, ntim, nlat*nlon);
[heat, clim] = remove_mean(heat);
default_global
biff = 3; num = 1;
cd /home/disk/hayes2/dvimont/csiro/matlab_data/Heat_Content
if biff == 1;
load RAW_L1-7_CEOF.mat; tit = 'Unfiltered Data';
elseif biff == 2;
load LP10_L1-7_CEOF.mat; tit = 'Lowpass Filtered Data ( > 10 Years )';
elseif biff == 3;
load HP8_L1-7_CEOF.mat; tit = 'Highpass Filtered Data ( < 8 Years )';
end
% Get time series
cwat = 4.218e3; % heat capacity of liquid water, J/(kg K)
rhowat = 1e3; % density of liquid water, kg/m^3
timeseries = sqrt(2)*pcs(:,num)./std(pcs(:,num));
lag = 0; lg = lag*pi/180;
num = 1; nfrm = 6;
cint = 50; clev = [-1000:cint:1000];
j = sqrt(-1);
figure(1); figure_orient;
for ind = 1:nfrm;
wgt = conj(exp(j * ((ind-1) * pi/nfrm + lg) ));
temtim = real(wgt .* timeseries);
tem = 1e-6*cwat*rhowat*temtim'*heat./ntim;
ccoef = corr_sig(temtim, heat);
tem = reshape(tem, nlat, nlon);
tstat = reshape((abs(ccoef) >= 0.25), nlat, nlon);
subplot(3,2,ind);
cla;
gshade2(abs(tstat), 0.9); hold on;
gcont(tem, clev); hold off;
dc;
title(['Phase = ' num2str((ind-1)*180/nfrm + lag)]);
if ind >= 5;
xlabel(['Cont. Int.: ' num2str(cint) 'x10^6 J m^-^2']);
end
end
subplot(3,2,3);
ylabel([tit ': 0:240m Heat Content Regressed on CPC1 of Heat Content; '...
'( r > 0.25 shaded )']);
cd /home/disk/tao/dvimont/matlab/CSIRO/Heat/Plot_Ceof
%print -dps2 HP8_HC_CPCHC_L1-7.ps
% Plot the CPC's 2222222222222222222222222222222222222222222222222222222222222
ind = 1;
if ind == 1; tim = 101:550; else; tim = 551:1000; end;
figure(2); figure_orient;
for m = 1:3
subplot(4,1,m)
tind = ((ntim/3)*(m-1)+[1:(ntim/3)]);
plot(tim(tind), real(timeseries(tind)), '-k', ...
tim(tind), imag(timeseries(tind)), '--k');
axis([min(tim(tind)-1) max(tim(tind)+1) -3 3]);
set(gca, 'YTick', [-3:3.5], 'XTick', [(tim(1)-1):10:max(tim)]);
grid
end
subplot(4,1,1)
title([tit ': CPC1 0:240m HEAT: ' num2str(round(per(1)))...
'% Variance Explained']);
subplot(4,1,2)
ylabel(['Solid = REAL, Dashed = IMAG']);
subplot(4,1,3)
xlabel('Years')
if 1;
subplot(4,2,7);
plot(real(timeseries), imag(timeseries), '.k');
axis square
axis([-4 4 -4 4]); grid on
xlabel('Re(CPC1)'); ylabel('Im(CPC1)');
end;
subplot(4,2,7);
binnum = 24;
binind = [(-1*(pi-(pi/binnum))):(2*pi/binnum):(pi-(pi/binnum))];
ph = atan2(imag(timeseries), real(timeseries));
hgram = hist(ph, binind);
bar((180/pi)*binind, hgram);
axis([-180 180 0 30])
set(gca, 'XTick', -180:90:180);
grid on
xlabel('PHASE (deg)');
ylabel('COUNT');
subplot(4,2,8); nind = 1;
tim = 1:450;
v1 = real(timeseries); v2 = imag(timeseries);
lab = ['Re(CPC1)'; 'Im(CPC1)'];
tind = (-8*nind):nind:(8*nind);
a = zeros(length(tind),1);
for lag = tind;
a((1/nind)*(lag+max(tind))+1) = corr(v1, v2, lag);
end
bar(tind, a);
title(['< ' lab(1,:) ', ' lab(2,:) ' >'])
% xlabel([lab(1,:) ' leads ' lab(2,:) ' | ' lab(2,:) ' leads ' lab(1,:)])
xlabel([lab(1,:) ' leads | ' lab(2,:) ' leads '])
axis([min(tind-1) max(tind+1) -1.1 1.1])
set(gca, 'XTick', min(tind):2*nind:max(tind), 'YTick', -1:.25:1)
grid on
%print -dps2 HP8_CPCHC_L1-7.ps