Documentation of ubar_grad_dHCprime

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

Help text

  Loop through HP LP and RAW

Cross-Reference Information

Listing of script ubar_grad_dHCprime

clear

for biff2 = 1:2

%  Start with ubardt'dx and vbardt'dy

biff = 7;
tim = 101:550;
ctlim = [106 304 -62.5 62.5];
lev = 1:biff;

cd /home/disk/hayes2/dvimont/csiro/data

%  First, get ubar and vbar

filin = 'ul_L1-10.nc';
nc = netcdf(filin, 'nowrite');
  latu = nc{'latitude'}(:);
  lonu = nc{'longitude'}(:);
  [xk, yk] = keep_var(ctlim, lonu, latu);
  ubar = nc{'u'}(tim,lev,yk,xk);
  mv = nc{'u'}.missing_value(:);
nc = close(nc);
ubar(ubar == mv) = NaN;
ubar = squeeze(mean2(ubar))/100;
latu = latu(yk); lonu = lonu(xk);

filin = 'vl_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;
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(ctlim, 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;

%  Load CPCs

cd /home/disk/hayes2/dvimont/csiro/matlab_data/Heat_Content
if biff2 == 1;
  load LP10_L1-7_CEOF.mat; tit = 'Lowpass Filtered Data ( > 10 Years )';
  ptit = 'LP10';
elseif biff2 == 2;
  load HP8_L1-7_CEOF.mat; tit = 'Highpass Filtered Data ( < 8 Years )';
  ptit = 'HP8';
elseif biff2 == 3;
  load RAW_L1-7_CEOF.mat; tit = 'Unfiltered Data';
  ptit = 'RAW';
end

%  Get regressions

lag = 0; lg = lag*pi/180; lg2 = 1;
num = 1; lind = 1;
nfrm = 6;

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

%  Plot the data

get_global; FRAME = [110 299 -60 60]; XAX = lonu; YAX = latu;
cwat = 4.218e3;  %  heat capacity of liquid water, J/(kg K)
rhowat = 1e3;    %  density of liquid water, kg/m^3

figure(1); figure_orient;
cint = 1; clev = [-50:cint:-cint cint:cint:50];
for i = 1:nfrm
  tem = cwat * rhowat * -1 * squeeze(udhdx(i,:,:));
  subplot(3,2,i);
    gcont(tem, clev);
    dc2(tem);
    title(['Phase = ' num2str((i-1)*180/nfrm + lag)]);
    box on
  if i > 4;
    xlabel(['Contour Interval:  ' num2str(cint) ' W m^-^2']);
  end
end
subplot(3,2,3);
  ylabel([tit ':  Zonal Advection of HC'' by Mean U;  '...
          'Depth = 0:' num2str(depth(biff)) 'm']);

cd /home/disk/tao/dvimont/matlab/CSIRO/Heat/Plot_HCadv
eval(['print -dps2 ' ptit '_ubar_dHCprimedx_yr1.ps']);

figure(2); figure_orient;
for i = 1:nfrm
  tem = cwat * rhowat * -1 * squeeze(vdhdy(i,:,:));
  subplot(3,2,i);
    gcont(tem, clev);
    dc2(tem)
    title(['Phase = ' num2str((i-1)*180/nfrm + lag)]);
    box on
  if i > 4;
    xlabel(['Contour Interval:  ' num2str(cint) ' W m^-^2']);
  end
end
subplot(3,2,3);
  ylabel([tit ':  Meridional Advection of HC'' by Mean V;  '...
          'Depth = 0:' num2str(depth(biff)) 'm']);

eval(['print -dps2 ' ptit '_vbar_dHCprimedy_yr1.ps'])

end