Documentation of get_wbar_dtprimedz

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

Function Synopsis

[wdhdz, lat_out, lon_out, depth_out] = ...

Help text

  [wdhdz, lat_out, lon_out, depth_out] = ...
          get_wbar_dtprimedz(pcs, lims, nfrm, tim, lev);

Cross-Reference Information

Listing of function get_wbar_dtprimedz

function [wdhdz, lat_out, lon_out, depth_out] = ...
          get_wbar_dtprimedz(pcs, lims, nfrm, tim, lev);

if nargin == 1;
  lims = [110 300 -60 60]; nfrm = 6; tim = 101:550; lev = 1:7;
elseif nargin == 2;
  nfrm = 6; tim = 101:550; lev = 1:7;
elseif nargin == 3;
  tim = 101:550; lev = 1:7;
elseif nargin == 4;
  lev = 1:7;
end

top_ocean = ismember(1, lev);

cdtem = ['cd ' eval('pwd')];

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

%  First, get wbar

filin = 'wl_L2-10.nc';
nc = netcdf(filin, 'nowrite');
  depthw = nc{'depth'}(:);
  latw = nc{'latitude'}(:);
  lonw = nc{'longitude'}(:);
  [xk, yk] = keep_var(lims, lonw, latw);
  if top_ocean
    wbar = nc{'wl'}(tim,lev,yk,xk);
  else
    wbar = nc{'wl'}(tim, [min(lev-1) lev], yk, xk);
  end
  mv = nc{'wl'}.missing_value(:);
nc = close(nc);
wbar(wbar == mv) = NaN;
wbar = squeeze(mean(wbar))/1e6;
latw = latw(yk); lonw = lonw(xk);
depthw = depthw/100;
if top_ocean
  we = getnc('we', lims, lev, tim);
else
  we = getnc('we', lims, [min(lev)-1 lev], tim);
end
we = squeeze(mean(we))/1e6;
wbar = wbar + we;

%  Now get tprime

cd /home/disk/hayes2/dvimont/csiro/data
filin = 'temp_L1-10.nc';
nc = netcdf(filin, 'nowrite');
  depth = nc{'depth'}(:);
  latt = nc{'latitude'}(:);
  lont = nc{'longitude'}(:);
  [xk, yk] = keep_var(lims, lont, latt);
  if top_ocean
    temp = nc{'temp'}(tim, [lev max(lev+1)], yk, xk);
  else
    temp = nc{'temp'}(tim, [min(lev)-[2 1] lev max(lev+1)], yk, xk);
  end
  mv = nc{'temp'}.missing_value(:);
nc = close(nc);
temp(temp == mv) = NaN;
latt = latt(yk); lont = lont(xk);
temp = squeeze(temp);

[temp, tbar] = remove_mean(temp);
[ntim, nlev, nlat, nlon] = size(temp);
depth = depth/100;

%  Assume PCS are sent in the command line

%cd /home/disk/hayes2/dvimont/csiro/matlab_data/Heat_Content
%load LP10_L1-7_CEOF.mat; tit = 'Lowpass Filtered Data ( > 10 Years )';
%load HP8_L1-7_CEOF.mat; tit = 'Highpass Filtered Data ( < 8 Years )';
%load RAW_L1-7_CEOF.mat; tit = 'Unfiltered Data';

%  Get regressions

if isreal(pcs);
  treg = regress_eof(temp, pcs, nfrm);
  lags = nfrm; nfrm = length(lags);
else
  treg = regress_ceof(temp, pcs, nfrm);
end

%  First, get dtdz
%  Interpolate to depthw, which is closer to temz.  This makes a smoother
%    interpolation, and will hopefully be more exact.

if top_ocean
  dz = diff(depth([lev max(lev+1)]));
  temz = depth(lev) + 0.5*dz;
else
  dz = diff(depth([min(lev)-[2 1] lev max(lev+1)]));
  temz = depth([min(lev)-[2 1] lev]) + 0.5*dz;
end
dtdz = reshape(shiftdim(treg, 1), nlev, nlat*nlon*nfrm);
dtdz = diff(dtdz) ./ (dz * ones(1, nlat*nlon*nfrm));
nlev = nlev - 1;
clear dtdz2
if top_ocean
  for i = 1:nlat*nlon*nfrm;
    dtdz2(:,i) = interp1([0; temz], [dtdz(1,i); dtdz(:,i)], depth(lev));
  end;
else
  for i = 1:nlat*nlon*nfrm;
    dtdz2(:,i) = interp1(temz, dtdz(:,i), depthw([min(lev-1) lev]));
  end;
end
nlev = size(dtdz2, 1);
dtdz = shiftdim(reshape(dtdz2, nlev, nlat, nlon, nfrm), 3);

%  Wait on this until after getting wdtdz
%  Interpolate wbar to new depths

%[nlev, nlat, nlon] = size(wbar);
%wbar = reshape(wbar, nlev, nlat*nlon);
%if top_ocean
%  for i = 1:nlat*nlon;
%    wbar2(:,i) = interp1([0; depthw(lev)], [0; wbar(:,i)], depth(lev));
%  end
%else
%  for i = 1:nlat*nlon;
%    wbar2(:,i) = interp1(depthw([min(lev-1) lev]), wbar(:,i), depth(lev));
%  end
%end
%nlev = length(lev);
%wbar = reshape(wbar2, nlev, nlat, nlon);

%  Calculate wbar_dtprimedz

wdtdz = (ones(nfrm, 1) * reshape(wbar, 1, nlev*nlat*nlon)) .* ...
         reshape(dtdz, nfrm, nlev*nlat*nlon);
wdtdz = reshape(wdtdz, nfrm, nlev, nlat, nlon);
wdtdz = reshape(shiftdim(wdtdz, 1), nlev, nlat*nlon*nfrm);

%  Interpolate to 'depth(lev)'

clear wdtdz2;
if top_ocean;
  for i = 1:nlat*nlon*nfrm;
    wdtdz2(:,i) = interp1([0; depthw(lev)], [0; wdtdz(:,i)], depth(lev));
  end
else
  for i = 1:nlat*nlon*nfrm;
    wdtdz2(:,i) = interp1(depthw([min(lev)-1 lev]), [wdtdz(:,i)], depth(lev));
  end
end
wdtdz = wdtdz2;
nlev = length(lev);

%  Get dz again

[tem1, tem2, depthw] = getll('wl', lims);
depthw = [0; depthw; 2*depth(10)-depthw(9)];
dz = diff(depthw([lev max(lev+1)]));

%  Integrate vertically

nlev = length(lev);
wdtdz = dz' * wdtdz;
wdhdz = shiftdim(reshape(wdtdz, nlat, nlon, nfrm), 2);

%  Convert units to heat

cwat = 4.218e3;  %  heat capacity of liquid water, J/(kg K)
rhowat = 1e3;    %  density of liquid water, kg/m^3

wdhdz = -1 * cwat * rhowat * wdhdz;
lat_out = latw;
lon_out = lonw;
depth_out = depth(lev);

eval(cdtem);