Documentation of w_intprime_dHCbardz

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

Help text

ctlim = [108 302 -62 62];

Cross-Reference Information

Listing of script w_intprime_dHCbardz

clear
for biff2 = 1:3;
tim = 101:550;
ctlim = [180 274 -7 7];
lev = 1:8;
biff = 7;

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

%  First get tbar

filin = 'temp_L1-10.nc';
nc = netcdf(filin, 'nowrite');
  depth = nc{'depth'}(:);
  latt = nc{'latitude'}(:);
  lont = nc{'longitude'}(:);
  [xk, yk] = keep_var(ctlim, lont, latt);
  temp = nc{'temp'}(tim, lev, yk, xk);
  mv = nc{'temp'}.missing_value(:);
nc = close(nc);
temp(temp == mv) = NaN;
latt = latt(yk); lont = lont(xk);
temp = squeeze(mean2(temp));

%  Now get u' and v'

filin = 'ul_L1-10.nc';
nc = netcdf(filin, 'nowrite');
  latu = nc{'latitude'}(:);
  lonu = nc{'longitude'}(:);
  [xk, yk] = keep_var(ctlim, lonu, latu);
  u = nc{'u'}(tim,lev,yk,xk);
  mv = nc{'u'}.missing_value(:);
nc = close(nc);
u(u == mv) = NaN;
latu = latu(yk); lonu = lonu(xk);
[u, climu] = remove_mean(u/100);

filin = 'vl_L1-10.nc';
nc = netcdf(filin, 'nowrite');
  v = nc{'v'}(tim,lev,yk,xk);
  mv = nc{'v'}.missing_value(:);
nc = close(nc);
v(v == mv) = NaN;
[v, climv] = remove_mean(v/100);

depth = depth/100;
u = squeeze(u); v = squeeze(v);
[ntim, nlev, nlat, nlon] = size(u);

%  Load CPCs

cd /home/disk/hayes2/dvimont/csiro/matlab_data/Heat_Content
if biff2 == 1;
  load LP10_L1-7_CEOF.mat; tit = 'Lowpass Filtered Data ( > 10 Years )';
  ptit = 'LP10';
elseif biff2 == 2;
  load HP8_L1-7_CEOF.mat; tit = 'Highpass Filtered Data ( < 8 Years )';
  ptit = 'HP8';
elseif biff2 == 3;
  load RAW_L1-7_CEOF.mat; tit = 'Unfiltered Data';
  ptit = 'RAW';
end

%  Get regressions

lag = 0; lg = lag*pi/180; lg2 = 1;
num = 1; lind = 1;
nfrm = 6;

%  Store all regressions under one variable

[nlev, nlat, nlon] = size(squeeze(climu));
j = sqrt(-1);
timeseries = sqrt(2)*pcs(:,num)./std(pcs(:,num));
u = reshape(u, ntim, nlev*nlat*nlon);
v = reshape(v, ntim, nlev*nlat*nlon);
clear temtim ureg vreg
for i = 1:nfrm
  wgt = conj(exp(j * ((i-1) * pi/(lg2*nfrm) + lg) ));
  temtim(:,i) = squeeze(real(wgt .* timeseries));
  ureg(i, :) = temtim(:,i)' * u ./ ntim;
  vreg(i, :) = temtim(:,i)' * v ./ ntim;
end
ureg = reshape(ureg, nfrm, nlev, nlat, nlon);
vreg = reshape(vreg, nfrm, nlev, nlat, nlon);

%  Get div(u)

clear divu lat2 lon2 j
for i = 1:nfrm;
  for j = 1:nlev;
    [divu(i,j,:,:), lat2, lon2] = sph_div1(squeeze(ureg(i,j,:,:)), ...
                                           squeeze(vreg(i,j,:,:)), latu, lonu, 0, 1);
  end
end

%  Get vvel (vvel @ lev x = sum(vvel(:,1:x,:,:)))
%  Assume variable 'depth' refers to the middle (vertically) of the layer to be integrated

[tem, nlev, nlat, nlon] = size(divu);
d1 = depth(1);
depthw = [];
for i = 1:nlev;
  depthw = [depthw; depth(i)+d1];
  d1 = depth(i+1) - depth(i) - d1;
end
dz = diff([0; depthw]);
dz = dz(1:nlev) * ones(1, nlat*nlon);
dz = reshape(dz, nlev, nlat, nlon);
clear vvel
for i = 1:nfrm;
  vvel(i,:,:,:) = squeeze(divu(i,:,:,:)) .* dz;
end

%  Interpolate to regular depth

depthw = [0; depthw(1:nlev)];
vvel = shiftdim(vvel, 1);
vvel = [zeros(size(vvel(1,:,:,:))); vvel];
vvel = reshape(vvel, nlev+1, nlat*nlon*nfrm);
vvel = cumsum(vvel);

clear vvel2
for i = 1:nlat*nlon*nfrm;
  vvel2(:,i) = interp1(depthw, vvel(:,i), depth(1:biff));
end
vvel2 = shiftdim(reshape(vvel2, biff, nlat, nlon, nfrm), 3);

%  Get dtdz

[nlev, nlat, nlon] = size(temp);
temp2 = reshape(temp, nlev, nlat*nlon);
dz = diff(depth(1:nlev)) * ones(1, nlat*nlon);
dzdepth = depth(1:(nlev-1)) + 0.5*diff(depth(1:nlev));
dtdz = diff(temp2) ./ dz;
dtdz = [dtdz(1,:); dtdz];
clear dtdz2;
for i = 1:nlat*nlon;
  dtdz2(:,i) = interp1([0; dzdepth], dtdz(:,i), depth(1:biff));
end

%  Integrate vertically

vvel2 = reshape(vvel2, nfrm, biff, nlat*nlon);
clear wdhdz
for i = 1:nfrm;
  wdhdz(i,:) = diff(depthw(1:(biff+1)))' * (squeeze(vvel2(i,:,:)) .* dtdz2);
end
wdhdz = reshape(wdhdz, nfrm, nlat, nlon);




%  Alternative Method

clear

tim = 101:550;
lims = [106 304 -62.5 62.5];
lev = 1:7;
nfrm = 6;
lag = 0;

% Load the PCs

for biff2 = 1:2;

cd /home/disk/hayes2/dvimont/csiro/matlab_data/Heat_Content
if biff2 == 1;
  load LP10_L1-7_CEOF.mat; tit = 'Lowpass Filtered Data ( > 10 Years )';
  ptit = 'LP10';
elseif biff2 == 2;
  load HP8_L1-7_CEOF.mat; tit = 'Highpass Filtered Data ( < 8 Years )';
  ptit = 'HP8';
elseif biff2 == 3;
  load RAW_L1-7_CEOF.mat; tit = 'Unfiltered Data';
  ptit = 'RAW';
end

%  Get wbdtdz

cd /home/disk/tao/dvimont/matlab/CSIRO/Heat
   [wdhdz, lat_out, lon_out, depth_out] = ...
           get_wprime_dtbardz(pcs, lims, nfrm, tim, lev);

%  Plot the data

default_global; XAX = lon_out; YAX = lat_out; FRAME = [110 299 -60 60];

figure(2); figure_orient;
cint = 1; clev = [-50:cint:-cint cint:cint:50];
for i = 1:nfrm
  tem = -1 * squeeze(wdhdz(i,:,:));
  subplot(3,2,i);
    gcont(tem, clev);
    dc2(tem)
    title(['Phase = ' num2str((i-1)*180/nfrm + lag)]);
  if i > 4;
    xlabel(['Contour Interval:  ' num2str(cint) ' W m^-^2']);
  end
end
subplot(3,2,3);
  ylabel([tit ':  Vertical Advection of HCbar by W'';  '...
          'Depth = 0:' num2str(depth_out(max(lev))) 'm']);

cd /home/disk/tao/dvimont/matlab/CSIRO/Heat/Plot_HCadv
eval(['print -dps2 ' ptit '_wprime_dHCbardz_yr1.ps']);

end