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_wprime_dtbardz(pcs, lims, nfrm, tim, lev);
| This function calls | |
|---|---|
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
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);
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(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);
wprime = nc{'wl'}(tim,lev,yk,xk);
mv = nc{'wl'}.missing_value(:);
nc = close(nc);
wprime(wprime == mv) = NaN;
latw = latw(yk); lonw = lonw(xk);
depthw = depthw/100;
wprime = wprime + getnc('we', lims, lev, tim);
[wprime, tem] = remove_mean(wprime/1e6);
[ntim, nlev, nlat, nlon] = size(wprime);
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
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));
wprime = reshape(wprime, ntim, nlev*nlat*nlon);
clear temtim wreg
for i = 1:nfrm
wgt = conj(exp(j * ((i-1) * pi/(lg2*nfrm) + lg) ));
temtim(:,i) = squeeze(real(wgt .* timeseries));
wreg(i, :) = temtim(:,i)' * wprime ./ ntim;
end
wreg = reshape(wreg, nfrm, nlev, nlat, nlon);
% First, get dtdz
[nlev, nlat, nlon] = size(temp);
dz = diff(depth([lev max(lev+1)]));
temz = depth(lev) + 0.5*dz;
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.
for i = 1:nlat*nlon;
dtdz2(:,i) = interp1([0; temz], [dtdz(1,i); dtdz(:,i)], depth(lev));
end;
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;
for i = 1:nlat*nlon*nfrm;
wdtdz2(:,i) = interp1([0; depthw(lev)], [0; wdtdz(:,i)], depth(lev));
end
wdtdz = wdtdz2;
% 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);