Documentation of sst_eof_heat_regress

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

Help text

 save sst_eof_bp_4.5-50yr.mat ld10 pc10 lam per kp lat lon ctlim

Cross-Reference Information

Listing of script sst_eof_heat_regress

clear
cd /home/disk/hayes2/dvimont/csiro/data
filin = ['temp_A_L1-10.nc'];
nc = netcdf(filin, 'nowrite');
  depth = nc{'depth'}(:);
  lat = nc{'latitude'}(:);
  lon = nc{'longitude'}(:);
  ctlim = [110 300 -75 65];
  [xk, yk] = keep_var(ctlim, lon, lat);
  temp = nc{'temp'}(:,1,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)));
lat = lat(yk); lon = lon(xk);
get_global; default_global; FRAME = ctlim;
temp = reshape(squeeze(temp), ntim, nlat*nlon);
temp = detrend(temp);
temp = cosweight(temp, lat);
clim = mean(temp);
kp = find(~isnan(clim)); npt = length(kp);
temp = temp(:, kp);
[b, a] = butter(6, (2/4.5));
[b2, a2] = butter(6, (2/50));
templp = filtfilt(b2, a2, temp);
tempbp = filtfilt(b, a, temp)-templp;
temphp = temp - filtfilt(b, a, temp);
c = tempbp*tempbp';
[lam, pcs, per] = eof(c);
pc10 = pcs(:,1:10);
wgt = std(pc10);
pc10 = (pc10 - ones(ntim, 1)*mean(pc10)) ./ (ones(ntim, 1)*wgt);
ld10 = pc10' * templp ./ ntim;
cd /home/disk/hayes2/dvimont/csiro/matlab_data

%  Now look at HEAT regressed on PC's of the above.

clear

cd /home/disk/hayes2/dvimont/csiro/data

filin = ['temp_A_L1-10.nc'];
nc = netcdf(filin, 'nowrite');

  depth = nc{'depth'}(:);
  lat = nc{'latitude'}(:);
  lon = nc{'longitude'}(:);

  ctlim = [110 300 -75 65];  %  Pacific only
%  ctlim = [-.1 360 -75 75];
  [xk, yk] = keep_var(ctlim, lon, lat);

  temp = nc{'temp'}(:,:,yk,xk);
  mv = nc{'temp'}.missing_value(:);

nc = close(nc);
[ntim, nlev, nlat, nlon] = size(temp);
clim = mean(squeeze(temp(:,7,:,:)));
tkp = find(reshape(clim, 1, nlat*nlon) ~= mv); nkp = length(tkp);
temp = reshape(temp, ntim, nlev, nlat*nlon);
temp = temp(:, :, tkp);
lat = lat(yk); lon = lon(xk);
get_global; default_global; FRAME = ctlim;

middepth = [0; (depth(1:7)+depth(2:8)) / 2];
wgt = diff(middepth)/100;

temp = shiftdim(temp, 1);
temp = reshape(temp, nlev, nkp*ntim);

heat = wgt' * temp(1:7,:);
heat = shiftdim(reshape(heat, nkp, ntim), 1);

temp = shiftdim(reshape(temp(1, :), nkp, ntim), 1);

temp = detrend(temp);
heat = detrend(heat);

%  Load top level EOF data

cd /home/disk/hayes2/dvimont/csiro/matlab_data

load sst_eof_lp_50yr.mat
pcLP = pc10; perLP = per;
%nlat = length(lat); nlon = length(lon);
%ntim = size(pc10, 1); npt = size(ld10, 1);

load sst_eof_bp_4.5-50yr.mat
pcBP = pc10; perBP = per;

load sst_eof_hp_4.5yr.mat
pcHP = pc10; perHP = per;

lab = ['LP'; 'BP'; 'HP'];
pol = -1*[1 1 -1];

c = 4.2*1000*27500*100*100*1e-12

num = 1;
for i = 1:3;

  if i == 1;
    tit = '50 Year Lowpass Filter Used for EOF';
  elseif i == 2;
    tit = '4.5 - 50 Year Bandpass Filter Used for EOF';
  elseif i == 3;
    tit = '4.5 Year Highpass Filter Used for EOF';
  end

  eval(['timtem = pc' lab(i,:) '(:,num) * pol(i);']);
  eval(['per = per' lab(i,:) '(num);']);

  tem1 = NaN * ones(1, nlat*nlon);
  tem2 = NaN * ones(1, nlat*nlon);
  tem1(tkp) = timtem' * temp ./ ntim; tem1 = reshape(tem1, nlat, nlon);
  tem2(tkp) = c*timtem' * heat ./ ntim; tem2 = reshape(tem2, nlat, nlon);
  tem1 = reshape(tem1, nlat, nlon);
  tem2 = reshape(tem2, nlat, nlon);

  figure(i);  figure_orient;

  sp(1);
    mcont(tem1, [-2:.05:2], 'giso');
    title(['<12.5m TEMP, PC' num2str(num) '\_' lab(i,:) '> : ' tit]);
    xlabel('Contour Interval:  0.05 C (std)^-^1');

  sp(2);
    mcont(tem2, [-200:10:200], 'giso');
    title(['<275m HEAT, PC' num2str(num) '\_' lab(i,:) '> : ' tit]);
    xlabel('Contour Interval: 10 * 10^1^2 J (std)^-^1');

  cd /home/disk/tao/dvimont/matlab/CSIRO/Plots
  eval(['print -dps2 TEMP12.5_HEAT_PC' num2str(num) '_' lab(i,:) '.ps']);

  figure(4);
  subplot(3,1,i);
    plot(1:1000, timtem);
    axis([0 1001 -3.5 3.5]);
    title(['PC' num2str(num) ' for ' lab(i,:) ' filtered 12.5m TEMP:  ' ...
           num2str(round(per)) '% Variance Explained']);
    set(gca, 'YTick', [-3:3]);
    xlabel('year');
    ylabel('std');
    grid on

end
  eval(['print -dps2 PC' num2str(num) '_LP_BP_HP.ps']);

cd /home/disk/tao/dvimont/matlab/CSIRO/Data
load butter_4.5_ctstar.mat