Global Index (short | long) | Local contents | Local Index (short | long)
First, find where the thermocline is
| This script calls | |
|---|---|
clear
tim = [101:550];
cd /home/disk/hayes2/dvimont/csiro/data
nc = netcdf('temp_A_L1-10.nc', 'nowrite');
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
depth = nc{'depth'}(:);
ctlim = [-.01 360 -90 90];
[xk, yk] = keep_var(ctlim, lon, lat);
temp = nc{'temp'}(tim, :, yk, xk);
mv = nc{'temp'}.missing_value(:);
nc = close(nc);
[ntim, nlev, nlat, nlon] = size(temp);
temp(find(temp == mv)) = NaN * ones(size(find(temp == mv)));
depth = depth / 100;
lat = lat(yk); lon = lon(xk);
% Find maximum of mean vertical temperature gradient.
if 0;
temp = reshape(temp, ntim, nlev*nlat*nlon);
[temp, clim] = remove_mean(temp);
clim = reshape(clim, nlev, nlat*nlon);
dtdz = diff(clim) ./ (-1 * diff(depth) * ones(1, nlat*nlon));
dtdz = rave(dtdz, 3);%,3),3);
depth2 = depth(1:nlev-1) + diff(depth)/2;
tcdepth = NaN * ones(1, nlat*nlon);
for i = 1:nlat*nlon;
if ~isnan(max(dtdz(:,i)));
ind = find(dtdz(:,i) == max(dtdz(:,i)));
if ind < 2; ind = 2; elseif ind > 8; ind = 8; end;
% tcdepth(i) = depth2(ind);
tcdepth(i) = ind;
end
end
% Calculate heat content above mean thermocline
temp = reshape(temp, ntim, nlev, nlat*nlon);
middepth = [0; (depth(1:9)+depth(2:10)) / 2];
wgt = diff(middepth);
temp2 = NaN*ones(ntim, nlat*nlon);
for i = 1:ntim;
for j = 1:nlat*nlon;
if ~isnan(tcdepth(j));
temp2(i,j) = squeeze(temp(i,(1:tcdepth(j)),j)) * wgt(1:tcdepth(j));
end
end;
end;
end
% Do the same thing as above, but for the instantaneous thermocline depth
temp = reshape(temp, ntim, nlev, nlat*nlon);
clim = squeeze(mean(mean(temp)));
kp = find(~isnan(clim)); nkp = length(kp);
temp = temp(:,:,kp);
temp = reshape(shiftdim(temp, 1), nlev, nkp*ntim);
depthdiff = -1 * diff(depth);
depth2 = depth(1:nlev-1) + diff(depth)/2;
middepth = [0; (depth(1:9)+depth(2:10)) / 2];
wgt = diff(middepth);
tcdepth = NaN * ones(1, nkp*ntim);
dtdz = rave((diff(temp)./ (depthdiff*ones(1,nkp*ntim))), 3);
[tem, tcdepth] = max(dtdz, [], 1);
tcdepth(find(tcdepth < 2)) = 2*ones(size(find(tcdepth < 2)));
tcdepth(find(tcdepth > 8)) = 8*ones(size(find(tcdepth > 8)));
heat = NaN * ones(1, nkp*ntim);
for i = 1:nkp*ntim;
heat(i) = wgt(1:tcdepth(i))' * temp((1:tcdepth(i)),i);
end
tc450 = tcdepth;
heat450 = heat;
cd /home/disk/hayes2/dvimont/csiro/matlab_data
%save heat_tc_yr101-550.mat heat450 tc450 kp lat lon depth depth2 tim
% Load new data
cd /home/disk/hayes2/dvimont/csiro/matlab_data
load heat_tc_yr101-550.mat
ntim = length(tim); nlat = length(lat); nlon = length(lon);
heat = NaN * ones(ntim, nlat*nlon);
heat(:,kp) = shiftdim(reshape(heat450, length(kp), ntim), 1);
heat = reshape(heat, ntim, nlat, nlon);
ctlim = [110 300 -30 30];
[xk, yk] = keep_var(ctlim, lon, lat);
lat = lat(yk); lon = lon(xk);
heat = heat(:,yk,xk);
tcdepth = NaN * ones(ntim, nlat*nlon);
tcdepth(:,kp) = shiftdim(reshape(tc450, length(kp), ntim), 1);
tcdepth = reshape(tcdepth, ntim, nlat, nlon);
tcdepth = tcdepth(:, yk, xk);
for i = 1:ntim;
for j = 1:length(yk);
for k = 1:length(xk);
if ~isnan(tcdepth(i,j,k));
tcdepth(i,j,k) = depth2(tcdepth(i,j,k));
end
end
end
end
% Try some stuff (like eofs) here
[ntim, nlat, nlon] = size(heat);
heat2 = NaN * ones(ntim, nlat/2, nlon/2);
lat2 = NaN * ones(1, nlat/2); lon2 = NaN * ones(1, nlon/2);
for i = 1:(nlat/2);
yind = 2*(i-1) + [1:2];
for j = 1:(nlon/2);
xind = 2*(j-1) + [1:2];
heat2(:,i,j) = squeeze(mean(mean(shiftdim(tcdepth(:,yind,xind),1))));
lon2(j) = mean(lon(xind));
end
lat2(i) = mean(lat(yind));
end
[ntim, nlat, nlon] = size(heat2);
heat2 = reshape(heat2, ntim, nlat*nlon);
clim = mean(heat2);
kp = find(~isnan(clim));
heat2 = detrend(heat2);
heat3 = heat2 ./ (ones(ntim, 1) * std(heat2));
heat3 = reshape(heat3, ntim, nlat, nlon);
heat3 = cosweight(heat3, lat2);
heat3 = reshape(heat3, ntim, nlat*nlon);
heat3 = heat3(:, kp);
[lam, per, lds, pcs] = eof_dan(heat3);
lat = lat2; lon = lon2;
get_global; default_global; FRAME = ctlim;
num = 1
tem = NaN * ones(1, nlat*nlon);
tem(kp) = lds(:,num);
tem = reshape(tem, nlat, nlon);
cint = .1; clev = [-1:cint:1];
figure(1)
sd(1)
gcont(tem, clev); dc;
sd(2)
plot(pcs(:,num));
tem = reshape(tcdepth, nlat, nlon);
figure(1); figure_orient; clf;
gcont(tem, depth2)
dc;
for i = 1:nlat
for j = 1:nlon
text(XAX(j), YAX(i), num2str(find(depth2 == tem(i,j))))
end
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots
figure(2); figure_orient;
for i = 2:8
ind = find(tcdepth == depth2(i));
% pat = mean2(shiftdim(clim(:,ind), 1))';
pat = clim(:,ind);
subplot(4,2,(i-1));
plot(pat, -1*depth, 'k'); hold on;
plot([5 30], -1*[depth2(i) depth2(i)], '-r'); hold off;
set(gca, 'YTick', [-550:100:-50], 'YTickLabel', [550:-100:50])
axis([5 30 -550 0])
grid on
title(['TC Depth = ' num2str(depth2(i))])
end