Documentation of gu_philander_thesis

Global Index (short | long) | Local contents | Local Index (short | long)

Help text

save gu_philand2.mat vbtp wbtp vptb wptb ubtp uptb lims tim lags

Cross-Reference Information

This script calls
clf corr dc2 default_global detrend fo frame gcont get_dof getct getheat getll getnc hline keep_var mean2 pncont regress_eof resid sptalk vert_ubtp_uptb vert_vdtdy_wdtdz vert_vptb_wptb vline

Listing of script gu_philander_thesis


clear
cd ~/matlab/CSIRO/Thesis/Data
load LP9_detrend_L1-7_EOF_yr101-1000.mat; pcs = 1*pcs;
lims = [178 227 -44 44];
tim = 101:1000;
lags = -15:15;
if 0;
cd ~/matlab/CSIRO/Thesis/Chap5
[vbtp, wbtp] = ...
              vert_vdtdy_wdtdz(pcs, lims, lags);
hello
[vptb, wptb] = ...
              vert_vptb_wptb(pcs, lims, lags);
hello
[ubtp, uptb] = vert_ubtp_uptb(pcs, lims, lags);
cd ~/matlab/CSIRO/Thesis/Data
end;


cd ~/matlab/CSIRO/Thesis/Data
load gu_philand.mat
%load gu_philand2.mat

temp = getnc('temp', lims, 1:10, [min(tim-1) tim]);
ntim = length(tim)+1;
dtdt = (temp(3:ntim,:,:,:) - temp(1:(ntim-2),:,:,:)) ./ (2*3600*24*365);
dtdt = detrend(dtdt);
[dtreg, dtcoef] = regress_eof(dtdt, -1*pcs(1:(ntim-2),1), lags);
[treg, tcoef] = regress_eof(temp(2:ntim,:,:,:), -1*pcs, lags);

[latt, lont, deptht] = getll('temp', lims);
[latv, lonv, depthv] = getll('v', lims);
[latw, lonw, depthw] = getll('wl', lims);
[latu, lonu, depthu] = getll('u', lims);

lims2 = [150 180 -44 44];
%lims2 = [180 225 -44 44];
[xku, yku] = keep_var(lims2, lonu, latu);
[xkw, ykw] = keep_var(lims2, lonw, latw);
[xkt, ykt] = keep_var(lims2, lont, latt);

ubtp = squeeze(mean2(shiftdim(ubtp(:,:,yku,xku), 3)));
vbtp = squeeze(mean2(shiftdim(vbtp(:,:,yku,xku), 3)));
wbtp = squeeze(mean2(shiftdim(wbtp(:,:,ykw,xkw), 3)));
uptb = squeeze(mean2(shiftdim(uptb(:,:,yku,xku), 3)));
vptb = squeeze(mean2(shiftdim(vptb(:,:,yku,xku), 3)));
wptb = squeeze(mean2(shiftdim(wptb(:,:,ykw,xkw), 3)));
dtreg = squeeze(mean2(shiftdim(dtreg(:,:,ykt,xkt), 3)));
treg = squeeze(mean2(shiftdim(treg(:,:,ykt,xkt), 3)));
dtcoef = squeeze(mean2(shiftdim(dtcoef(:,:,ykt,xkt), 3)));
tcoef = squeeze(mean2(shiftdim(tcoef(:,:,ykt,xkt), 3)));

[nlag, nlev, nlat, nlon] = size(vbtp);

for i=1:nlag;
  ubtp2(i,:,:) = interp2(latu, depthu, squeeze(ubtp(i,:,:)), latw', depthw);
  vbtp2(i,:,:) = interp2(latv, depthv, squeeze(vbtp(i,:,:)), latw', depthw);
  uptb2(i,:,:) = interp2(latu, depthu, squeeze(uptb(i,:,:)), latw', depthw);
  vptb2(i,:,:) = interp2(latv, depthv, squeeze(vptb(i,:,:)), latw', depthw);
  treg2(i,:,:) = interp2(latt, deptht, squeeze(dtreg(i,:,:)), latw', depthw);
end

uptb3 = uptb2 + vptb2 + wptb;
ubtp3 = ubtp2 + vbtp2 + wbtp;

tem1 = 1e9*(treg2-uptb3-ubtp3); tit1 = 'Resid';
tem2 = 1e9*ubtp3; tit2 = 'UBTP';
lat1 = latw; depth1 = depthw;

biff = 5;
tem1 = 1e9*ubtp3; tit1 = ['\bfU\rm\nabla\cdotT'''];
tem2 = 1e9*uptb3; tit2 = ['\bfu\rm''\nabla\cdotT'];
lat1 = latw; depth1 = depthw;

tem1 = 4*dtcoef; tit1 = 'T''_t';
tem2 = 4*tcoef; tit2 = 'T''';
lat1 = latt; depth1 = deptht;

tem2 = 1e9*(treg2 - uptb3 - ubtp3);





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lagplot = -8:1:-4; lp = length(lagplot);

for biff = 1:5;

if biff == 1;
tem2 = 1e9*dtreg; tit2 = 'T''_t'; cint2 = 10; 
  unit2 = 'x 10^-^8 K s^-^1 std^-^1';
tem1 = 20*treg; tit1 = 'T'''; cint1 = 0.05; 
  unit1 = 'K std^-^1';
lat1 = latt; depth1 = deptht;
elseif biff == 2;
tem1 = 1e9*ubtp; tit1 = 'UT''_x'; cint1 = 10; 
  unit1 = 'x 10^-^8 K s^-^1 std^-^1';
tem2 = 1e9*uptb; tit2 = 'u''T_x'; cint2 = 10; 
  unit2 = 'x 10^-^8 K s^-^1 std^-^1';
lat1 = latu; depth1 = depthu;
elseif biff == 3;
tem1 = 1e9*vbtp; tit1 = 'VT''_y'; cint1 = 10; 
  unit1 = 'x 10^-^8 K s^-^1 std^-^1';
tem2 = 1e9*vptb; tit2 = 'v''T_y'; cint2 = 10; 
  unit2 = 'x 10^-^8 K s^-^1 std^-^1';
lat1 = latu; depth1 = depthu;
elseif biff == 4;
tem1 = 1e9*wbtp; tit1 = 'WT''_z'; cint1 = 10; 
  unit1 = 'x 10^-^8 K s^-^1 std^-^1';
tem2 = 1e9*wptb; tit2 = 'w''T_z'; cint2 = 10; 
  unit2 = 'x 10^-^8 K s^-^1 std^-^1';
lat1 = latw; depth1 = depthw;
elseif biff == 5;
tem1 = 1e9*ubtp3; tit1 = ['\bfU\rm\nabla\cdotT''']; cint1 = 10; 
  unit1 = 'x 10^-^8 K s^-^1 std^-^1';
tem2 = 1e9*uptb3; tit2 = ['\bfu\rm''\nabla\cdotT']; cint2 = 10; 
  unit2 = 'x 10^-^8 K s^-^1 std^-^1';
lat1 = latw; depth1 = depthw;
end

%  Plot the data;

figure(biff); fo(1); clf;

for i = 1:lp;
if lagplot(i) < -4;
  cint = 0.2; clev = sort([-cint:-cint:min(min(min([tem1; tem2]))) ...
                           cint:cint:max(max(max([tem1; tem2])))]);
else
  cint = 0.2; clev = sort([-cint:-cint:min(min(min([tem1; tem2]))) ...
                           cint:cint:max(max(max([tem1; tem2])))]);
end
  sptalk(6, 2, 2*i-1);
    tind = find(lags == lagplot(i));
    tem = squeeze(1*tem1(tind,:,:)); 
    pncont(lat1, -1*depth1, tem, clev, 0, 'k');
    set(gca, 'YTick', -500:125:0, 'YTickLabel', 500:-125:0, ...
             'XTick', -40:10:40);
    axis([-40 40 -545 0]);
    yl(i) = ylabel(['Lag = ' num2str(lagplot(i))]);
    hold on;
      vline(-40:10:40, ':k');
      hline(-500:125:-125, ':k');
    hold off;
    if i ~= lp;
%      set(gca, 'XTickLabel', []);
    end
    if i == 1;
      ti(1) = title(tit1);
    elseif i == lp;
      xl(1) = xlabel(['Contour Interval:  ' num2str(cint*cint1) ' ' unit1]);
    end
    set(gca, 'fontsize', 9, 'box', 'on');
end

for i = 1:lp;
if lagplot(i) < -4;
  cint = 0.2; clev = sort([-cint:-cint:min(min(min([tem1; tem2]))) ...
                           cint:cint:max(max(max([tem1; tem2])))]);
else
  cint = 0.2; clev = sort([-cint:-cint:min(min(min([tem1; tem2]))) ...
                           cint:cint:max(max(max([tem1; tem2])))]);
end
  sptalk(6, 2, 2*i);
    tind = find(lags == lagplot(i));
    tem = squeeze(1*tem2(tind,:,:)); 
    pncont(lat1, -1*depth1, tem, clev, 0, 'k');
    set(gca, 'YTick', -500:125:0, 'YTickLabel', 500:-125:0, ...
             'XTick', -40:10:40);
    axis([-40 40 -545 0]);
    hold on;
      vline(-40:10:40, ':k');
      hline(-500:125:-125, ':k');
    hold off;
    if i ~= lp;
%      set(gca, 'XTickLabel', []);
    end
%    set(gca, 'YTickLabel', []);
    if i == 1;
      ti(2) = title(tit2);
    elseif i == lp;
      xl(2) = xlabel(['Contour Interval:  ' num2str(cint*cint2) ' ' unit2]);
    end
    set(gca, 'fontsize', 9, 'box', 'on');
end

set(yl, 'fontsize', 11);
set(xl, 'fontsize', 10);
set(ti, 'fontsize', 10);

end;  %  biff loop





cd ~/matlab/CSIRO/Thesis/Chap5/Plots
if 0;
  figure(1); print -dps2 GP_ttend_temp_neglag.ps
  figure(2); print -dps2 GP_uterms_neglag.ps
  figure(3); print -dps2 GP_vterms_neglag.ps
  figure(4); print -dps2 GP_wterms_neglag.ps
  figure(5); print -dps2 GP_ubtp_uptb_neglag.ps
end






%%%%%%%%%%%%%%%  Look at lagged autocorrelation

clear

ct = getct;
np = getnc('temp', [180 205 25 37.5], 1, 101:1000);
%np = getheat([150 180 5 20], 5:8, 101:1000);
np = squeeze(mean2(mean2(shiftdim(np, 1))));

[b, a] = butter(9, 2/9);
ct = filtfilt(b, a, detrend(ct));
np = filtfilt(b, a, detrend(np));

tind = [1:450];
tind = [451:900];
tind = [1:900];

clear c ctc
lags = -15:15;
for i = 1:length(lags);
  c(i) = corr(np(tind), ct(tind), lags(i));
  ctc(i) = corr(ct(tind), ct(tind), lags(i));
end

figure(1); clf;
sptalk(4,2,1);
  h = bar(lags, c);
  axis([-13 13 -0.6 0.4])
  set(gca, 'XTick', -10:5:10, 'YTick', -1:.2:1);
  grid on
  set(h, 'facecolor', [.7 .7 .7]);
  title('< NP\_SST, CT >');
  xlabel('NP\_SST leads CT   |   CT leads NP\_SST')

if 0;
lind = find(lags == 0);
sptalk(4,2,2);
  h = bar(lags, ctc.*c(lind));
  axis([-13 13 -0.6 0.4])
  set(gca, 'XTick', -10:5:10, 'YTick', -1:.2:1);
  grid on
  set(h, 'facecolor', [.7 .7 .7]);
  title('CT autocorrelation'); 
end

%  Look at fishz stat:  Get expected correlation


rho = ctc*c(lind);
sigz = 1/sqrt(0.5*(get_dof(ct)+get_dof(np))-3);
z = 0.5*log((1+c)./(1-c));
muz = 0.5*log((1+rho)./(1-rho));

sptalk(4,2,2)
  h = bar(lags, z);
  hold on
  plot(lags, muz, '-k', lags, muz+2.0*sigz, '--k', ...
       lags, muz-2.0*sigz, '--k');
  hold off
  axis([-13 13 -0.9 0.5])
  set(gca, 'XTick', -10:5:10, 'YTick', -1:.2:1);
  grid on
  set(h, 'facecolor', [.7 .7 .7]);
  title('Fisher''s Z-statistic for NP\_SST and CT');
  xlabel('NP\_SST leads CT   |   CT leads NP\_SST')

cd ~/Thesis/Chap5
%print -dps2 GP_fishz_stat.ps




%%%%%%%%%%%%%%%  Look at regression on NP

clear

np = getnc('temp', [180 205 25 37.5], 1, 101:1000);
np = mean2(mean2(shiftdim(np, 1)));

lims2 = [150 180 -44 44];
%lims2 = [180 225 -44 44];
tim = 101:1000;

temp = getnc('temp', lims2, 1:10, tim);

temp = detrend(temp);
[b, a] = butter(9, 2/9);
temp = filtfilt(b, a, temp);
np = filtfilt(b, a, np);
temp = squeeze(mean2(shiftdim(temp, 3)));

lags = -6:16;
[reg1, c1] = regress_eof(temp, np', lags);
[lat, lon, depth, lm] = getll('temp', lims2);

lagplot = [-6:2:16];
for i = 1:length(lagplot);
lind = find(lags == lagplot(i));
sptalk(6,2,i);
pncont(lat, -depth, squeeze(reg1(lind,:,:)), [-.5:0.025:.5], 0, 'k');
axis([-40 40 -500 0])
ylabel([num2str(lagplot(i))]);
end

hc = getheat([100 300 -45 45], 4:7, 101:1000);
hc = detrend(hc);
hc = filtfilt(b, a, hc);

[reg2, c2] = regress_eof(hc, np', lags);
[lat, lon, depth, lm] = getll('temp', [100 300 -45 45]);

figure(2); fo(1); clf;
default_global; FRAME = [105 295 -40 40];

lagplot = [-6:2:16];
for i = 1:length(lagplot);
lind = find(lags == lagplot(i));
sptalk(6,2,i);
gcont(1e-8*reg2(i,:,:), 0.25);
dc2(lm);
ylabel([num2str(lagplot(i))]);
end