Global Index (short | long) | Local contents | Local Index (short | long)
Get regressions of fluxes:
| This script calls | |
|---|---|
clear
cd /home/disk/hayes2/dvimont/csiro/matlab_data
load LPregmap_3deof_-67.5to60pac.mat
ctlim = FRAME;
get_global; default_global; FRAME = ctlim;
ntim = size(pcs, 1);
[nlev, nlat, nlon] = size(clim);
tim = [101:550];
cd /home/disk/hayes2/dvimont/csiro/data
filin = ['temp_A_L1-10.nc'];
nc = netcdf(filin, 'nowrite');
lind = 1:10;
depth = nc{'depth'}(lind);
nc = close(nc);
cd /home/disk/hayes2/dvimont/csiro/data
nc = netcdf('zonal_stress_A_1000_years.nc', 'nowrite');
lat1 = nc{'latitude'}(:);
lon1 = nc{'longitude'}(:);
[xk, yk] = keep_var(FRAME, lon1, lat1);
taux = nc{'smfzon'}(tim,yk,xk);
mv = nc{'smfzon'}.missing_value(:);
nc = close(nc);
taux(find(taux == mv)) = NaN * ones(size(find(taux == mv)));
nc = netcdf('merid_stress_A_1000_years.nc', 'nowrite');
tauy = nc{'smfmer'}(tim,yk,xk);
mv = nc{'smfmer'}.missing_value(:);
nc = close(nc);
tauy(find(tauy == mv)) = NaN * ones(size(find(tauy == mv)));
lat1 = lat1(yk); lon1 = lon1(xk);
nc = netcdf('heat_flux_A_1000_years.nc', 'nowrite');
lat2 = nc{'latitude'}(:);
lon2 = nc{'longitude'}(:);
[xk, yk] = keep_var(FRAME, lon2, lat2);
hflx = nc{'stfht'}(tim,yk,xk);
mv = nc{'stfht'}.missing_value(:);
nc = close(nc);
hflx(find(hflx == mv)) = NaN * ones(size(find(hflx == mv)));
lat2 = lat2(yk); lon2 = lon2(xk);
nc = netcdf('temp_A_L1-10.nc', 'nowrite');
lat = nc{'latitude'}(:);
lon = nc{'longitude'}(:);
[xk, yk] = keep_var(FRAME, lon, lat);
for i = 1:4;
temp(:,i,:,:) = nc{'temp'}(tim,(3*(i-1)+1),yk,xk);
end
mv = nc{'temp'}.missing_value(:);
nc = close(nc);
temp(find(temp == mv)) = NaN * ones(size(find(temp == mv)));
lat = lat(yk); lon = lon(xk);
depth = depth(1:3:10);
hflx2 = hflx;
hflx = hflx2 + (5.67e-8 * (squeeze(temp(:,1,:,:)+273.15).^4));
hflx = reshape(hflx, ntim, size(hflx,2)*size(hflx,3));
taux = reshape(taux, ntim, size(taux,2)*size(taux,3));
tauy = reshape(tauy, ntim, size(tauy,2)*size(tauy,3));
[hflx, clim] = remove_mean(hflx);
[taux, clim] = remove_mean(taux); kpx = find(~isnan(clim));
[tauy, clim] = remove_mean(tauy); kpy = find(~isnan(clim));
%taux = taux(:,kpx); tauy = tauy(:,kpy);
[b, a] = butter(6, 2/4.5);
hflx = filtfilt(b, a, hflx);
taux = filtfilt(b, a, taux);
tauy = filtfilt(b, a, tauy);
hflxr = pcs.' * hflx ./ ntim;
tauxr = pcs.' * taux ./ ntim;
tauyr = pcs.' * tauy ./ ntim;
hflxr = reshape(hflxr, 10, size(hflx, 2), size(hflx, 3));
tauxr = reshape(tauxr, 10, size(taux, 2), size(taux, 3));
tauyr = reshape(tauyr, 10, size(tauy, 2), size(tauy, 3));
% Plot these in the following format:
% Phase0 Phase45 Phase90 Phase135
% level1
% level4
% level7
% level10
% taux,tauy
% hflx
cint1 = [0.1 0.1 0.05 0.025 1.2 1.5];
for i = 1:6;
fig1 = floor(i/2 - 1/2);
if i < 5;
[ntim, nlev, nlat, nlon] = size(temp);
tem1 = reshape(squeeze(regmap(1,(3*(i-1)+1),:,:)), 1, nlat*nlon);
tem2 = std(reshape(squeeze(temp(:,i,:,:)), ntim, nlat*nlon));
cint = cint1(i); clev = [-1:cint:1];
xlab = [num2str(cint) ' K (std)^-^1'];
XAX = lon; YAX = lat;
tit = ['Depth: ' num2str(depth(i)/100) ' m; '];
elseif i == 5;
tem1 = tauxr(1,:);
tem2 = tauyr(1,:);
cint = cint1(i);
XAX = lon1; YAX = lat1;
nlat = length(YAX); nlon = length(XAX);
tit = ['Wind Stress, '];
elseif i == 6;
[ntim, nlat, nlon] = size(hflx);
tem1 = -1*hflxr(1,:);
tem2 = std(hflx);
cint = cint1(i); clev = [-18:cint:18];
xlab = [num2str(cint) ' W m^-^2 (std)^-^1'];
XAX = lon2; YAX = lat2;
nlat = length(YAX); nlon = length(XAX);
tit = ['ATM Heat Flux (Pos. Up), '];
end;
dofx = ntim/4.5;
if i ~= 5; % Perform regular t-test on correllation coefficient
lag = 0; lg = lag*pi/180;
num = 1; nfrm = 6;
for ind = 1:nfrm;
fig2 = floor(ind/2 - 0.5);
wgt = conj(exp(j * ((ind-1) * pi/nfrm + lg) ));
tem = squeeze(real(wgt .* tem1));
ccoef = tem./tem2;
tstat = ccoef ./ (sqrt((1-ccoef.^2)/dofx));
score = tscore(dofx, 2.5);
tem = reshape(tem, nlat, nlon);
tstat = reshape(tstat, nlat, nlon);
figure(3*fig1 + fig2+1); figure_orient;
subplot(2,2,(rem(ind-1,2)+1+2*rem(i-1,2)));
cla;
gshade2(abs(tstat), score); hold on;
gcont(tem, clev); hold off;
dc;
title([tit 'Phase = ' num2str((ind-1)*180/nfrm + lag)]);
xlabel(['Cont. Int.: ' xlab]);
end
else
lag = 0; lg = lag*pi/180;
num = 1; nfrm = 6;
for ind = 1:nfrm;
fig2 = floor(ind/2 - 0.5);
wgt = conj(exp(j * ((ind-1) * pi/nfrm + lg) ));
tem3 = reshape(squeeze(real(wgt .* tem1)), nlat, nlon);
tem4 = reshape(squeeze(real(wgt .* tem2)), nlat, nlon);
figure(3*fig1 + fig2+1); figure_orient;
subplot(2,2,(rem(ind-1,2)+1+2*rem(i-1,2)));
% figure(2*fig1 + fig2+1); figure_orient;
% subplot(2,2,(rem(ind-1,2)+1+2*rem(i-1,2)));
[h, a, c] = gquiv(tem3, tem4, 1.2); dc;
title([tit 'Phase = ' num2str((ind-1)*180/nfrm + lag)]);
xlabel(['Max = ' num2str(round(c(2)*100)/100) ' dynes cm^-^1'])
end
end
end
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots4
% Plot time series
figure(1); figure_orient;
for m = 1:4
subplot(5,1,m)
tind = ((ntim/4)*(m-1)+[1:(ntim/4)]);
plot((100+tind), real(pcs(tind, num)), '-k', ...
(100+tind), imag(pcs(tind, num)), '--k');
axis([100+min(tind-1) 100+max(tind+1) -3 3]);
set(gca, 'YTick', [-3:3], 'XTick', [100:10:(ntim+100)]);
grid
end
subplot(5,1,1)
title(['LP\_PAC 3D TEMP (-67.5 to 60): CPC', num2str(num), ...
': Dashed = IMAG, Solid = REAL']);
subplot(5,1,4)
xlabel('Years')
tim = 1:450;
v1 = real(pcs(tim,1)); v2 = imag(pcs(tim,1));
lab = ['Re(CPC1)'; 'Im(CPC1)'];
tind = -10:10;
a = zeros(length(tind),1);
for lag = tind;
a(lag+max(tind)+1) = corr(v1, v2, lag);
end
subplot(5,2,9);
bar(tind, a);
title(['LP: < ' 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])
set(gca, 'XTick', min(tind):2:max(tind), 'YTick', -1:.25:1)
grid on
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots4
% Plot the phase of the time series vs. amplitude
ph = atan2(imag(pcs(:,1)), real(pcs(:,1)));
amp = sqrt(imag(pcs(:,1)).^2 + real(pcs(:,1)).^2);
wgt = (4 / (2*pi));
figure(1); figure_orient;
for m = 1:4
subplot(5,1,m)
tind = ((ntim/4)*(m-1)+[1:(ntim/4)]);
plot((100+tind), amp(tind), '-k', ...
(100+tind), wgt.*(ph(tind)+pi), '.k');
axis([100+min(tind-1) 100+max(tind+1) 0 4]);
set(gca, 'YTick', [0:4], 'XTick', [100:10:(ntim+100)]);
set(gca, 'YTickLabel', [-180:90:180]);
grid
end
subplot(5,1,1)
title(['LP\_PAC 3D TEMP (-67.5 to 60): CPC', num2str(num), ...
': Solid = AMPLITUDE, Dots = PHASE']);
subplot(5,1,4)
xlabel('Years')
subplot(5,2,9);
plot(real(pcs(:,1)), imag(pcs(:,1)), '.');
axis square; axis([-4 4 -4 4]);
grid on
xlabel('REAL(CPC1)'); ylabel('IMAG(CPC1)');
subplot(5,2,10);
binnum = 24;
binind = [(-1*(pi-(pi/binnum))):(2*pi/binnum):(pi-(pi/binnum))];
hgram = hist(ph, binind);
bar((180/pi)*binind, hgram);
axis([-180 180 0 35])
xlabel('PHASE (deg)');
ylabel('COUNT');
cd /home/disk/tao/dvimont/matlab/CSIRO/Plots4