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