% calculate reference pressure and altitude
% assumes netCDF file containing MM5 output has already been opened
% code based on getrefheight.m/getrefpressure.m by Matt Wyant
%
% G Mullendore
%
% z = z at the sigma levels (full)
% zh = z at thermodynamic point (half)
% p0, p0h = same for pressure

if ~exist('header','var')
%get ptop, psfc, Tsfc, lapse
  load_MM5header
end

R = 287.04;
g = 9.81;

sigmah = fid{'sigma_level'}(:);
nzh = length(sigmah);
sigmah_mat = repmat(sigmah,[1 ny-1 nx-1]);
sigma = fid{'sigma_level_full'}(:);
sigma_mat = repmat(sigma,[1 ny-1 nx-1]);
nz = length(sigma);

fprintf('This code assumes that stratospheric temperature is zero K\n');
  %see getrefheight.m for details
pstarcrs = fid{'pstarcrs'}(:); %(REFERENCE) SURFACE PRESSURE MINUS PTOP
pstarcrs_mat = shiftdim(repmat(pstarcrs,[1 1 nz]),2);
pstarcrsh_mat = shiftdim(repmat(pstarcrs,[1 1 nz-1]),2);

%p0 = (psfc-ptop)*sigma+ptop;
%p0h = (psfc-ptop)*sigmah+ptop;
p0 = pstarcrs_mat .* sigma_mat +ptop;
p0h = pstarcrsh_mat .* sigmah_mat +ptop;

clear pstar*

%compute T0, T0h
T0 = zeros(size(p0));
T0(end,:,:)=Tsfc;
T0h = zeros(size(p0h));
for k = nzh:-1:1,
  T0(k,:,:) = T0(k+1,:,:) - lapse*log(p0(k+1,:,:)./p0(k,:,:));
  T0h(k,:,:) = T0(k+1,:,:) - lapse*log(p0(k+1,:,:)./p0h(k,:,:));
end

%compute z starting at the surface and climbing upwards using
%hydrostatic relation
z = zeros(size(p0));
terrain = fid{'terrain'}(:);
z(end,:,:) = terrain;
for k = nz-1:-1:1,
  z(k,:,:) = z(k+1,:,:) - (R/g) * T0h(k,:,:) .* log(p0(k,:,:)./p0(k+1,:,:));
end
zh = zeros(size(p0h));
zh = .5 * (z(1:end-1,:,:)+z(2:end,:,:));

