Documentation of get_wprime_dtbardz

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_wprime_dtbardz(pcs, lims, nfrm, tim, lev);

Cross-Reference Information

Listing of function get_wprime_dtbardz

function [wdhdz, lat_out, lon_out, depth_out] = ...
          get_wprime_dtbardz(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

%  Now get tbar

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(mean(temp));

%  Next, get wprime

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
    wprime = nc{'wl'}(tim,lev,yk,xk);
  else
    wprime = nc{'wl'}(tim, [min(lev-1) lev], yk, xk);
  end
  mv = nc{'wl'}.missing_value(:);
nc = close(nc);
wprime(wprime == mv) = NaN;
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
wprime = wprime + we;

[wprime, tem] = remove_mean(wprime/1e6);
[ntim, nlev, nlat, nlon] = size(wprime);
clear we;

wprime = reshape(wprime, ntim, nlev*nlat*nlon);

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);
  wreg = regress_eof(wprime, pcs, nfrm);
  lags = nfrm; nfrm = length(lags);
else
  wreg = regress_ceof(wprime, pcs, nfrm);
end

%  First, get dtdz

[nlev, nlat, nlon] = size(temp);
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(temp, nlev, nlat*nlon);
dtdz = diff(dtdz) ./ (dz * ones(1, nlat*nlon));
nlev = nlev - 1;

%  As this isn't a very smooth function, interpolate to depthw (which is 
%    very near temz, so the interpolation should be close).  Interpolate
%    wdhdz later.

clear dtdz2
if top_ocean
  for i = 1:nlat*nlon;
    dtdz2(:,i) = interp1([0; temz], [dtdz(1,i); dtdz(:,i)], depth(lev));
  end;
else
  for i = 1:nlat*nlon;
    dtdz2(:,i) = interp1(temz, dtdz(:,i), depthw([min(lev-1) lev]));
  end;
end
nlev = size(dtdz2, 1);
dtdz = reshape(dtdz2, nlev, nlat, nlon);

%  Wait on this:
%  Interpolate wreg to new depths

%wreg = reshape(shiftdim(wreg, 1), nlev, nlat*nlon*nfrm);
%clear wreg2;
%for i = 1:nlat*nlon*nfrm;
%  wreg2(:,i) = interp1([0; depthw(lev)], [0; wreg(:,i)], depth(lev));
%end
%wreg = shiftdim(reshape(wreg2, nlev, nlat, nlon, nfrm), 3);

%  Calculate wprime_dtbardz

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

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;

%  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

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);