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
| This script calls | |
|---|---|
clear
filin = 'temp_A_L1-10.nc';
% filin = 'temp_M_L5_1000_years_new.nc'
cd /home/disk/hayes2/dvimont/csiro/data
tim = 101:550;%
%tim = round(12*[240:1/12:(279+11/12)]+1);
if 0;
nc = netcdf(filin, 'nowrite');
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
ctlim = [110 300 -60 60];
% ctlim = [115 285 -30 30];
[xk, yk] = keep_var(ctlim, lon, lat);
temp = nc{'temp'}(tim,1,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);
else
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);
end
if 1;
[ntim, nlat, nlon] = size(temp);
temp = shiftdim(temp, 1);
for i = 1:nlat/2
indy = 2*(i-1)+[1:2];
for j = 1:nlon/2
indx = 2*(j-1)+[1:2];
temp2(:,i,j) = squeeze(mean(mean(temp(indy,indx,:))));
end
end
lat = mean([lat(1:2:nlat)'; lat(2:2:nlat)']);
lon = mean([lon(1:2:nlon)'; lon(2:2:nlon)']);
temp = temp2; clear temp2;
[ntim, nlat, nlon] = size(temp);
end;
temp = reshape(temp, ntim, nlat*nlon);
%temp = rave(temp, 5);
temp = reshape(temp, ntim, nlat, nlon);
get_global;
XAX = lon; YAX = lat; FRAME = ctlim;
%[temp, clim] = annave(temp);
%temp2 = (temp - ones(ntim, 1)*mean2(temp));
temp = reshape(temp, ntim, nlat*nlon);
[temp, clim] = remove_mean(temp);
%clim = mean(temp);
temp = detrend(temp);
temp = cosweight(temp, lat);
kp = find(~isnan(clim));
temp = reshape(temp, ntim, nlat*nlon);
temp = temp(:, kp);
% Low Pass Filter the data:
[b,a]=butter(6,2/4.5);
temp = temp - filtfilt(b, a, temp);
%timekp = 11:1:(ntim-10);
%temp = temp(timekp, :);
[ntim, nkp] = size(temp);
temp = detrend(temp);
[lam, lds, pcs, per] = complex_eof(temp);
cd /home/disk/tao/dvimont/matlab/CSIRO/Data
% save ceof_T1_yr101-550.mat lam lds pcs per kp lat lon
if 0;
tem = lds;
lds = pcs;
pcs = tem;
end
pc10 = pcs;
lds10 = lds;
j = sqrt(-1);
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(nkp, 1)*weight1) + ...
j*imag(lds).*(ones(nkp, 1)*weight2);
cd /home/disk/tao/dvimont/matlab/CSIRO/Data
% save ceof_L1-10.mat pcs lds kp lat lon nlat nlon ntim
load ceof_L1-10.mat
figure(1)
num = 2;
tem = NaN * ones(1, nlat*nlon);
tem(kp) = real(lds(:,num))';
tem = reshape(tem, nlat, nlon);
sp(1);
if ismap(gca); clma; end;
gcont(tem, [-1:.05:1]);%, 'giso');
dc
title(['AVERAGE TEMP to 620m: REAL CEOF' num2str(num)])
xlabel(['Contour Interval: 0.025 C (std)^-^1'])
tem = NaN * ones(1, nlat*nlon);
tem(kp) = imag(lds(:,num))';
tem = reshape(tem, nlat, nlon);
sp(2);
if ismap(gca); clma; end;
gcont(tem, [-1:.025:1]);%, 'giso');
dc
title(['AVERAGE TEMP to 620m: IMAG CEOF' num2str(num)])
xlabel(['Contour Interval: 0.025 C (std)^-^1'])
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots
figure(2); orient tall
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(['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
% Take a look at phase space
figure(3)
ind = 1:ntim;
axis([-3 3 -3 3])
x = real(pcs(ind, num)); y = imag(pcs(ind, num));
x = x - rave(x, 50); y = y - rave(y, 50);
x = interp(x, 50);
y = interp(y, 50);
cometnew(x, y, .05);
plot(x, y, '.');
ctstar = squeeze(mean2(mean2(shiftdim(temp, 1))));
clear ctann
%ctann(1) = mean(ctstar(1:2));
for i = 1:((ntim/12));
ctann(i) = mean(ctstar(12*(i-1) + [1:12]));
end
nx = length(ctann);
figure(2)
for i = 1:2
if i==1;
ctann = real(pcs(:,1))';
elseif i==2;
ctann = imag(pcs(:,1))';
end
nx = length(ctann);
% Remove trend from ctann:
%x = 1:nx; x = (x-mean(x)) / std(x);
%a1 = ctann * x' / nx;
%ctann = ctann - a1*x;
%%ctann = ctann - rave(ctann, 10);
%ctann = ctann(2:5:1000);
% Look at spectral characteristics of ctstar:
nfft = 32;
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)) + ...
% (corr(ctann, ctann, 3)^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));
rner = rn * 1.67;
dofx = nx / n2;
sp(i);
semilogy(f2, p(:,1), 'b-', f2, rn, 'b-', f2, rner, 'b--')
axis([0 0.5 0.1 15])
grid
end
ctlow = detrend(clow);
% Look at phases of lds:
cd /home/disk/hayes2/dvimont/csiro/data
tim = 101:550;
nc = netcdf(filin, 'nowrite');
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
ctlim = [110 300 -60 60];
[xk, yk] = keep_var(ctlim, lon, lat);
temp = nc{'temp'}(tim,1,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);
temp2 = temp; XAX = lon1; YAX = lat1;
nlat = length(YAX); nlon = length(XAX);
temp2 = detrend(temp2);
FRAME = [ctlim(1:2) floor(min(YAX)) ceil(max(YAX))];
lag = 0;
figure(1); clf; figure_orient;
lg = lag*pi/180; lg2 = 1;
num = 1; lind = 1;
nfrm = 6;
cint = 0.1; clev = [-9:cint:9];
for i = 1:nfrm;
wgt = conj(exp(j * ((i-1) * pi/(lg2*nfrm) + lg) ));
temtim = squeeze(real(wgt .* pcs(:,num)));
temtim = (temtim - mean(temtim)) ./ std(temtim);
tem = temtim' * temp ./ ntim;
tem = reshape(tem, nlat, nlon);
subplot(3,2,i);
gcont(tem, clev); hold off; dc;
title(['Phase = ' num2str(((i-1)*180)/(lg2*nfrm)+lag)]);
xlabel(['Contour Interval: ' num2str(cint) ' K std^-^1'])
end
% Plot time series
temtim = sqrt(2)*pcs(:,num)./std(pcs(:,num));
figure(2); figure_orient;
for m = 1:3
subplot(4,1,m)
tind = ((ntim/3)*(m-1)+[1:(ntim/3)]);
plot(tim(tind), real(temtim(tind)), '-k', ...
tim(tind), imag(temtim(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 20 DEG ISOTHERM: ' num2str(round(per(1)))...
'% Variance Explained']);
subplot(4,1,2)
ylabel(['Solid = REAL, Dashed = IMAG']);
subplot(4,1,3)
xlabel('Years')
if 0;
subplot(4,2,7);
plot(real(temtim), imag(temtim), '.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(pcs(:,1)), real(pcs(:,1)));
hgram = hist(ph, binind);
bar((180/pi)*binind, hgram);
axis([-180 180 0 40])
set(gca, 'XTick', -180:90:180);
grid on
xlabel('PHASE (deg)');
ylabel('COUNT');
subplot(4,2,8); nind = 1;
tim = 1:450;
v1 = real(pcs(tim,1)); v2 = imag(pcs(tim,1));
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
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots4