Documentation of get_ubar_gradtprime

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[udhdx, vdhdy, lat_out, lon_out, depth_out] = ...

Help text

  [ub_dhpdx, vb_dhpdy, lat_out, lon_out, depth_out] = ...
          get_ubar_gradtprime(pcs, lims, nfrm, tim, lev);

Cross-Reference Information

Listing of function get_ubar_gradtprime

function [udhdx, vdhdy, lat_out, lon_out, depth_out] = ...
          get_ubar_gradtprime(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 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

if isreal(pcs);
  treg = regress_eof(temp, pcs, nfrm);
  lags = nfrm; nfrm = length(lags);
else
  treg = regress_ceof(temp, pcs, nfrm);
end

%  Get thickness of layer, to convert t' to Heat'

[tem1, tem2, depthw] = getll('wl', lims);
depthw = [0; depthw; 2*depth(10)-depthw(9)];
dz = diff(depthw([lev max(lev+1)]));
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);