Global Index (short | long) | Local contents | Local Index (short | long)
[wdhdz, lat_out, lon_out, depth_out] = ...
[wdhdz, lat_out, lon_out, depth_out] = ...
get_wbar_dtprimedz(pcs, lims, nfrm, tim, lev);
| This function calls | |
|---|---|
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
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);
wbar = nc{'wl'}(tim,lev,yk,xk);
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;
we = getnc('we', lims, lev, tim);
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);
temp = nc{'temp'}(tim, [lev max(lev+1)], yk, xk);
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
lag = 0; lg = lag*pi/180; lg2 = 1;
num = 1; lind = 1;
% Store all regressions under one variable
j = sqrt(-1);
timeseries = sqrt(2)*pcs(:,num)./std(pcs(:,num));
temp = reshape(temp, ntim, nlev*nlat*nlon);
clear temtim treg
for i = 1:nfrm
wgt = conj(exp(j * ((i-1) * pi/(lg2*nfrm) + lg) ));
temtim(:,i) = squeeze(real(wgt .* timeseries));
treg(i, :) = temtim(:,i)' * temp ./ ntim;
end
treg = reshape(treg, nfrm, nlev, nlat, nlon);
% First, get dtdz
% Interpolate to depthw, which is closer to temz. This makes a smoother
% interpolation, and will hopefully be more exact.
dz = diff(depth([lev max(lev+1)]));
temz = depth(lev) + 0.5*dz;
dtdz = reshape(shiftdim(treg, 1), nlev, nlat*nlon*nfrm);
dtdz = diff(dtdz) ./ (dz * ones(1, nlat*nlon*nfrm));
nlev = nlev - 1;
for i = 1:nlat*nlon*nfrm;
dtdz2(:,i) = interp1([0; temz], [dtdz(1,i); dtdz(:,i)], depthw(lev));
end;
dtdz = shiftdim(reshape(dtdz2, nlev, nlat, nlon, nfrm), 3);
% Wait on this until after getting wdtdz
% Interpolate wbar to new depths
%wbar = reshape(wbar, nlev, nlat*nlon);
%for i = 1:nlat*nlon;
% wbar2(:,i) = interp1([0; depthw(lev)], [0; wbar(:,i)], depth(lev));
%end
%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;
for i = 1:nlat*nlon*nfrm;
wdtdz2(:,i) = interp1([0; depthw(lev)], [0; wdtdz(:,i)], depth(lev));
end
% Get dz again
d1 = depth(1);
depth2 = [];
for i = 1:nlev;
depth2 = [depth2; depth(i)+d1];
d1 = depth(i+1) - depth(i) - d1;
end
dz = diff([0; depth2]);
% 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);