Global Index (short | long) | Local contents | Local Index (short | long)
The following files are in the directory below: filin = 'temp_A_L1_1000_years.nc' % SST filin = 'temp_A_L1-5_1000_years.nc' % Averaged Pot. Temp., Top 160m filin = 'temp_A_L1-10_1000_years.nc' % Averaged Pot. Temp., Top 620m filin = 'temp_M_L1_1000_years_new.nc' % Monthly SST filin = 'temp_M_L5_1000_years_new.nc' % Monthly layer 5 temp
| This script calls | |
|---|---|
clear
cd /home/disk/hayes2/dvimont/csiro/data
nc = netcdf(filin, 'nowrite');
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
% ctlim = [130 285 -50 65];
ctlim = [180 270 -6 6];
% ctlim = [-.1 360 -90 90];
[xk, yk] = keep_var(ctlim, lon, lat);
temp = nc{'temp'}(:,1,yk,xk);
% temp = nc{'temp'}(:,yk,xk);
mv = nc{'temp'}.missing_value(:);
nc = close(nc);
temp = squeeze(temp);
[ntim, nlat, nlon] = size(temp);
temp(find(temp == mv)) = NaN * ones(size(find(temp == mv)));
lat = lat(yk); lon = lon(xk);
%temp = reshape(temp, ntim, nlat*nlon);
%temp = rave(temp, 1);
%temp = reshape(temp, ntim, nlat, nlon);
get_global;
XAX = lon; YAX = lat; FRAME = ctlim;
if 0
ctlim2 = [180 270 -6 6];
[xk, yk] = keep_var(ctlim2, lon, lat);
ctstar = squeeze(mean2(mean2(shiftdim(temp(:,yk,xk),1))));
end
%ctann = detrend(ctstar);
%ctann = (ctann - rave(rave(ctann, 2), 2));
%ctann = rave(rave(ctann, 2), 2);
ctann = detrend(chi);
ctann = (ctann - mean(ctann)) / std(ctann);
ctpat = ctann' * reshape(temp, ntim, nlat*nlon) / ntim;
tem = reshape(ctpat, nlat, nlon);
%[temp, clim] = annave(temp);
%temp2 = (temp - ones(ntim, 1)*mean2(temp));
%[temp, clim] = remove_mean(temp);
temp = reshape(temp, ntim, nlat*nlon);
clim = mean2(temp);
temp = detrend(temp);
temp = reshape(temp, ntim, nlat, nlon);
temp = cosweight(temp, lat);
kp = find(~isnan(clim));
temp = reshape(temp, ntim, nlat*nlon);
temp = temp(:, kp);
c = temp' * temp / ntim;
[lam, lds, per] = eof(c);
lds10 = lds(:,1:10);
lds10 = (lds10 - ones(size(kp'))*mean(lds10)) ./ (ones(size(kp'))*std(lds10));
pc10 = temp * lds10 / length(kp);
pc10 = pc10 ./ (ones(ntim, 1) * std(pc10));
lds10 = pc10' * temp / ntim;
figure(1)
num = 1;
tem = NaN * ones(1, nlat*nlon);
tem(kp) = lds10(num,:);
tem = reshape(tem, nlat, nlon);
sd(1);
gshade(tem, [-1:.05:1]);
colorbar
dc
sd(2);
plot(pc10(:,num));
% plot(ctann);
set(gca, 'XTick', [0:100:1000])
axis([0 1001 -3 3])
grid
ctstar = squeeze(mean2(mean2(shiftdim(temp, 1))));
[ctstar, clim] = annave(ctstar);
if 0
clear ctann
ctann(1) = mean(ctstar(1:2));
for i = 2:((ntim/3));
ctann(i) = mean(ctstar(3*(i-1) + [0:2]));
end
ctann = detrend(ctann);
nx = length(ctann);
end
ctann = detrend(pc10(:,1)');
nx = length(ctann);
if 1;
% ctann = detrend(ctstar(101:550));
ctann = detrend(ctstar(551:1000));
% ctann = detrend(ctann);
nx = length(ctann);
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Data
load ceof_heat1-7_yr551-1000lp.mat
load ceof_heat1-7_yr101-550lp.mat
ctannlp = real(pcslp(:,1)); fig = 3;
ctannlp = imag(pcslp(:,1)); fig = 4;
load ceof_heat1-7_yr551-1000.mat
load ceof_heat1-7_yr101-550.mat
ctann = real(pcslp(:,1));
ctann = imag(pcslp(:,1));
% Look at spectral characteristics of ctstar:
nfft = 128;
noverlap = nfft/2;
[p, f] = spectrum(ctannlp, nfft, noverlap);
n2 = nfft/2;
f2 = 0.5 * (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(fig)
sd(1)
semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', ...
f2, rner5, 'b--', f2, rner1, 'b-.')
% axis([0 0.5 0.005 1])
% axis([0 0.5 0.05 10])
% axis([0 2 .001 5])
axis([0 0.2 1e4 1e7])
grid
title('Power Spectrum for CSIRO CT: 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');
sd(2)
plot(1:length(ctann), ctann);
title('Detrended CT index, Annual average');
xlabel('Years')
ylabel('Degrees C')
grid
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots
% Butterworth Filter:
yr = [100 200];
subplot(3,2,6);
for i = 1:2
[b,a]=butter(6,(2/(yr(i))));
[h,w]=freqz(b,a,1024);
r=abs(h);
rlow=r.^2;
f2=w/(2*pi);
i
if i == 2;
rlow2 = rlow;
else
rlow1 = rlow;
end
end;
plot(f2, rlow1, '--k', f2, rlow2, '-k');
set(gca, 'box', 'on');
grid on
axis([0 .5 -.1 1.1])
title('Butterworth Filter response: N=6')
legend('4 year', '4.5 year');
clow = filtfilt(b, a, ctstar);
chi = ctstar - clow;
ctann = clow;
ctann = chi;
% Plot the data
ctann = detrend(clow);
for j = 1:2;
if j == 1;
ctann = ctstar; holdon = 0;
elseif j == 2;
ctann = ct12; holdon = 1;
end
for i = 1:4;
subplot(4,1,i);
if holdon; hold on; end;
ind = (250*(i-1)) + [1:250];
a = plot(ind, ctann(ind), '-k');
if holdon; set(a, 'linewidth', 1); end;
grid on
axis([250*(i-1) 250*i -.7 .7])
set(gca, 'YTick', [-1:.5:1]);
if holdon; hold off; end;
end
end
xlabel('Year');
subplot(4,1,1)
title('Lowpass filtered CT Index (4.5 year cutoff)')
subplot(3,2,5);
set(gca, 'XTick', [0:.1:.5])
subplot(3,2,6);
xlabel('Frequency: yr^-^1')
% Define indices
ctlow = rave(rave(detrend(ctstar), 3), 3);
ctlow = detrend(clow);
ctlow = (ctlow - mean(ctlow)) / std(ctlow);
ctpat = ctlow' * reshape(temp, ntim, nlat*nlon) / ntim;
ctlowpat = reshape(ctpat, nlat, nlon);
cthi = detrend(chi);
cthi = (cthi - mean(cthi)) / std(cthi);
ctpat = cthi' * reshape(temp, ntim, nlat*nlon) / ntim;
cthipat = reshape(ctpat, nlat, nlon);
ctdetrend = detrend(ctstar);
ctdetrend = (ctdetrend - mean(ctdetrend)) / std(ctdetrend);
ctpat = ctdetrend' * reshape(temp, ntim, nlat*nlon) / ntim;
ctdetrendpat = reshape(ctpat, nlat, nlon);
%cd /home/disk/tao/dvimont/matlab/CSIRO/Data
%save ct_csiro.mat ctlow cthi ctdetrend ctlowpat cthipat ctdetrendpat
% Plot some regression patterns:
ctlow = cthi;
ctlowpat = cthipat;
cint = .05;
clims = sort(unique([0:-cint:min(min(ctlowpat)) 0:cint:max(max(ctlowpat))]));
clims = [min(clims)-cint clims max(clims)+cint];
figure(2)
sd(1)
gshade(ctlowpat, clims)
colorbar
title('CSIRO SST regressed on Lowpass Filtered CT Index');
xlabel('Contour Interval: 0.05 C (std)^-^1');
dc
sd(2)
plot(1:1000, ctlow)
title('Detrended, Low Pass Filtered CT index, Annual Data');
xlabel('Years')
ylabel('STD')
grid
axis([0 1001 -3 3])
set(gca, 'YTick', [-3:3]);
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots