Documentation of look_at_stoc_run


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


Help text

cd /home/disk/hayes2/dvimont/ocean/coup/coup

Cross-Reference Information

This script calls

Listing of script look_at_stoc_run


clear
global RADUS
ctlim = [80 160 -6 6];
cd /home/disk/hayes2/dvimont/ocean/coup/stoc
nc = netcdf('sal.cdf', 'nowrite');
  lon = nc{'x'}(:);
  lat = nc{'y'}(:);

  lat = 1000*lat*180/(pi*RADUS);
  lon = 1000*lon*180/(pi*RADUS);

  [xk, yk] = keep_var(ctlim, lon, lat);
  cth = nc{'h',1}(:,:,yk,xk);
  cth = squeeze(cth);
  
nc = close(nc);

cth = squeeze(mean(mean(shiftdim(cth, 1))))-15000;

figure(2); fo(1); clf
subplot(4,1,1);
  hh = plot([1:1200]./12, -rave(cth, 3), 'k', ...
            [1:1200]./12, -rave(rave(cth, 25), 37), '-k');
  set(hh(2), 'linewidth', 2);
  set(gca, 'XTick', [1:120:1200]/12, 'XTickLabel', [0:10:100]);
  axis([0 100 -1250 1250]);
  set(gca, 'YTick', -5000:500:5000)
  grid on;
  
subplot(4,2,3);
  tem = cth(1:600);
  nfft = 256; nover = nfft/2; n2 = nfft/2;
  [p,f] = spectrum(tem, nfft, nover);
  rho = (corr(tem, tem, 1)+sqrt(corr(tem, tem, 2))+...
	 corr(tem, tem, 3)^(1/3))/3;
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn(3:(n2+1)));
  rn = rn * mean(p(3:(n2+1),1));
  rner1 = rn * 2.5; rner5 = rn * 1.92;
  dofx = length(tem) / n2;
  xt = [128 64 32 16 8 4 2 1 0.5 0.25];
  f2 = f.*6;
  hh = plot(log(f2), f2.*p(:,1), '-k', ...
	    log(f2), f2.*rn', '-k', ...
	    log(f2), f2.*rner5', '--k', ...
	    log(f2), f2.*rner1', '-.k');
%  hh = loglog(f2, p(:,1), '-k', f2, rn, '-k', ...
%	      f2, rner5, '--k', f2, rner1, '-.k');
  set(hh(1), 'linewidth', 2);
  set(gca, 'XTick', log(1./xt), 'XTickLabel', xt);
  axis([log(1/20) log(6.5) 1e6*[.1 5]]);
%  axis([1/20 6.1 5e2 5e6]); set(gca, 'YTick', 10.^[1:8]);
  grid on
  title('Years 1-50');
  
subplot(4,2,5);
  tem = cth(601:1200);
  nfft = 256; nover = nfft/2; n2 = nfft/2;
  [p,f] = spectrum(tem, nfft, nover);
  rho = (corr(tem, tem, 1)+sqrt(corr(tem, tem, 2))+...
	 corr(tem, tem, 3)^(1/3))/3;
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn(3:(n2+1)));
  rn = rn * mean(p(3:(n2+1),1));
  rner1 = rn * 2.5; rner5 = rn * 1.92;
  dofx = length(tem) / n2;
  xt = [128 64 32 16 8 4 2 1 0.5 0.25];
  f2 = f.*6;
  hh = plot(log(f2), f2.*p(:,1), '-k', ...
	    log(f2), f2.*rn', '-k', ...
	    log(f2), f2.*rner5', '--k', ...
	    log(f2), f2.*rner1', '-.k');
  set(hh(1), 'linewidth', 2);
  set(gca, 'XTick', log(1./xt), 'XTickLabel', xt);
  axis([log(1/20) log(6.5) 1e6*[-0.25 2.25]]);
  grid on
  title('Years 51-100');
  
subplot(4,2,7);
  tem = cth(1:1200);
  nfft = 256; nover = nfft/2; n2 = nfft/2;
  [p,f] = spectrum(tem, nfft, nover);
  rho = (corr(tem, tem, 1)+sqrt(corr(tem, tem, 2))+...
	 corr(tem, tem, 3)^(1/3))/3;
  rn = (1-rho^2) ./ (1-2*rho*cos((pi/n2)*(0:n2))+rho^2);
  rn = rn / mean(rn(3:(n2+1)));
  rn = rn * mean(p(3:(n2+1),1));
  rner1 = rn * 1.7; rner5 = rn * 2.5;
  dofx = length(tem) / n2;
  xt = [128 64 32 16 8 4 2 1 0.5 0.25];
  f2 = f.*6;
  hh = plot(log(f2), f2.*p(:,1), '-k', ...
	    log(f2), f2.*rn', '-k', ...
	    log(f2), f2.*rner5', '--k', ...
	    log(f2), f2.*rner1', '-.k');
  set(hh(1), 'linewidth', 2);
  set(gca, 'XTick', log(1./xt), 'XTickLabel', xt);
  axis([log(1/20) log(6.5) 1e6*[-0.25 2.25]]);
  grid on
  title('Years 1-100');
  

%  Plot regressions and such

nc = netcdf('sal.cdf', 'nowrite');
  h = nc{'h',1}(:,:,:,:);
nc = close(nc);
h = squeeze(h-15000);
h = h(:,:,2:3:180);

h = detrend(h);

[ntim, nlat, nlon] = size(h);
h = reshape(h, ntim, nlat*nlon);
h2 = rave(rave(h, 25), 37);

h = reshape(h, ntim, nlat, nlon);
h2 = reshape(h2, ntim, nlat, nlon);

cth2 = rave(rave(cth, 25), 37);


cd /home/disk/hayes2/dvimont/ocean/coup/coup
%save h_lph_thinned_stoc.mat h h2 cth cth2 lat lon

load h_lph_thinned_stoc.mat

[reg1, c1] = regress_eof(h, cth, 1);
[reg2, c2] = regress_eof(h-h2, cth-cth2, 1);
[reg3, c3] = regress_eof(h2, cth2, 1);

default_global; FRAME = [0 160 -45 45]; XAX = lon(2:3:180);
subplot(4,2,4); cla;
  gcont(reg1(1,:,:), 250, 1);
  title('Unfiltered')
  ylabel('2.5m Contours')
  set(gca, 'XTick', [0:40:160], 'XTickLabel', [-80:40:80]);

subplot(4,2,6); cla;
  gcont(reg2(1,:,:), 250, 1);
  title('5yr HP')
  ylabel('2.5m Contours')
  set(gca, 'XTick', [0:40:160], 'XTickLabel', [-80:40:80]);
  
subplot(4,2,8); cla;
  gcont(reg3(1,:,:), 100, 1);
  title('5yr LP')
  ylabel('1m Contours')
  set(gca, 'XTick', [0:40:160], 'XTickLabel', [-80:40:80]);

  
cd ~/matlab/Ocean/Figs
print -dps2 stoc_regs_ct.ps

  
figure(2); fo(1); clf;
nlag = length(-3:3);
for i = 1:4;
  subplot(4,2,(2*i)-1);
    gcont(reg2(i,:,:), 200, 1);
  subplot(4,2,(2*i));
    gcont(reg3(i,:,:), 50, 1);
end
figure(3); fo(1); clf;
nlag = length(-3:3);
for i = 1:3;
  subplot(4,2,(2*i)-1);
    gcont(reg2(i+4,:,:), 200, 1);
  subplot(4,2,(2*i));
    gcont(reg3(i+4,:,:), 25, 1);
end


figure(2); fo(1); clf;
[xk1, yk1] = keep_var([-1 161 -20 -15]);
tem = squeeze(mean(h(:,yk1,xk1), 2));
subplot(1,3,1);
  surface(lon(2:3:180), [1:1200]./12, tem);
  shading flat
  axis([0 160 0 100])
  set(gca, 'XTick', 0:40:160)
  caxis([-3000 3000])
[xk1, yk1] = keep_var([-1 161 -10 -5]);
tem = squeeze(mean(h(:,yk1,xk1), 2));
subplot(1,3,2);
  surface(lon(2:3:180), [1:1200]./12, tem);
  shading flat
  axis([0 160 0 100])
  set(gca, 'XTick', 0:40:160)
  caxis([-3000 3000])
[xk1, yk1] = keep_var([-1 161 -5 5]);
tem = squeeze(mean(h(:,yk1,xk1), 2));
subplot(1,3,3);
  surface(lon(2:3:180), [1:1200]./12, tem);
  shading flat
  axis([0 160 0 100])
  set(gca, 'XTick', 0:40:160)
  caxis([-3000 3000])

figure(3); fo(1); clf;
[xk1, yk1] = keep_var([-1 161 -20 -15]);
tem = squeeze(mean(h2(:,yk1,xk1), 2));
subplot(1,3,1);
  surface(lon(2:3:180), [1:1200]./12, tem);
  shading flat
  axis([0 160 0 100])
  set(gca, 'XTick', 0:40:160)
  caxis([-1000 1000])
[xk1, yk1] = keep_var([-1 161 -10 -5]);
tem = squeeze(mean(h2(:,yk1,xk1), 2));
subplot(1,3,2);
  surface(lon(2:3:180), [1:1200]./12, tem);
  shading flat
  axis([0 160 0 100])
  set(gca, 'XTick', 0:40:160)
  caxis([-1000 1000])
[xk1, yk1] = keep_var([-1 161 -5 5]);
tem = squeeze(mean(h2(:,yk1,xk1), 2));
subplot(1,3,3);
  surface(lon(2:3:180), [1:1200]./12, tem);
  shading flat
  axis([0 160 0 100])
  set(gca, 'XTick', 0:40:160)
  caxis([-1000 1000])