Global Index (short | long) | Local contents | Local Index (short | long)
[udhdx, vdhdy, lat_out, lon_out, depth_out] = ...
[ub_dhpdx, vb_dhpdy, lat_out, lon_out, depth_out] = ... get_ubar_gradtprime(pcs, lims, nfrm, tim, lev);
This function calls | |
---|---|
function [udhdx, vdhdy, lat_out, lon_out, depth_out] = ... get_ubar_gradtprime(pcs, lims, lags, tim, lev); nfrm = length(lags); 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 ubar and vbar filin = 'u_L1-10.nc'; nc = netcdf(filin, 'nowrite'); latu = nc{'latitude'}(:); lonu = nc{'longitude'}(:); [xk, yk] = keep_var(lims, lonu, latu); ubar = nc{'u'}(tim,lev,yk,xk); mv = nc{'u'}.missing_value(:); nc = close(nc); ubar(ubar == mv) = NaN; filin = 'ue_L1-10.nc'; nc = netcdf(filin, 'nowrite'); yk2 = [yk; max(yk)+1]; ue = nc{'ue'}(tim,lev,yk2,xk); mv = nc{'ue'}.missing_value(:); nc = close(nc); ue(ue == mv) = NaN; [ntim, nlev, nlat, nlon] = size(ue); ue = (ue(:,:,1:(nlat-1),:) + ue(:,:,2:nlat,:))/2; ubar = ubar + ue; clear ue; ubar = squeeze(mean2(ubar))/100; latu = latu(yk); lonu = lonu(xk); filin = 'v_L1-10.nc'; nc = netcdf(filin, 'nowrite'); vbar = nc{'v'}(tim,lev,yk,xk); mv = nc{'v'}.missing_value(:); nc = close(nc); vbar(vbar == mv) = NaN; filin = 've_L1-10.nc'; nc = netcdf(filin, 'nowrite'); xk2 = [xk; max(xk)+1]; ve = nc{'ve'}(tim,lev,yk,xk2); mv = nc{'ve'}.missing_value(:); nc = close(nc); ve(ve == mv) = NaN; [ntim, nlev, nlat, nlon] = size(ve); ve = (ve(:,:,:,1:(nlon-1)) + ve(:,:,:,2:nlon))/2; vbar = vbar + ve; clear ve; vbar = squeeze(mean2(vbar))/100; % 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, 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 treg = regress_eof(temp, pcs, lags); %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); % Get thickness of layer, to convert t' to Heat' d1 = depth(1); depthw = []; for i = 1:nlev; depthw = [depthw; depth(i)+d1]; d1 = depth(i+1) - depth(i) - d1; end dz = diff([0; depthw]); dz = shiftdim(reshape((dz * ones(1, nlat*nlon*nfrm)), nlev, nlat, nlon, nfrm), 3); treg = treg.*dz; % Get dtdx and dtdy global DEGREE RADIAN RADUS clear ty tx dtdy dtdx j for i = 1:nfrm; for lind = 1:nlev; for j = 1:nlat; ty(j,:) = interp1(lont, squeeze(treg(i,lind,j,:)), lonu)'; end for j = 1:nlon; tx(:,j) = interp1(latt, squeeze(treg(i,lind,:,j)), latu); end [dtdy(i,lind,:,:) temlat] = sph_grady1(ty, RADIAN*latt, RADIAN*lonu, 0); [dtdx(i,lind,:,:) temlon] = sph_gradx1(tx, RADIAN*latu, RADIAN*lont, 0); end end % Get ubar*dtdx and vbar*dtdy, and sum vertically [nlev, nlat, nlon] = size(ubar); ubar = reshape(ubar, 1, nlev*nlat*nlon); dtdx = reshape(dtdx, nfrm, nlev*nlat*nlon); udhdx = dtdx .* (ones(nfrm, 1) * ubar); udhdx = shiftdim(reshape(udhdx, nfrm, nlev, nlat, nlon), 1); udhdx = shiftdim(squeeze(sum(udhdx)), 2); vbar = reshape(vbar, 1, nlev*nlat*nlon); dtdy = reshape(dtdy, nfrm, nlev*nlat*nlon); vdhdy = dtdy .* (ones(nfrm, 1) * vbar); vdhdy = shiftdim(reshape(vdhdy, nfrm, nlev, nlat, nlon), 1); vdhdy = shiftdim(squeeze(sum(vdhdy)), 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 vdhdy = cwat * rhowat * vdhdy; udhdx = cwat * rhowat * udhdx; lat_out = latu; lon_out = lonu; depth_out = depth(lev); eval(cdtem);