%function [x] = getlev(aray,lev);
%This function takes a 3-dimensional array (aray) with values depending on 
%sigma level, latitude, and longitude, and interpolates the values to
%one pressure level (lev), generating a 2-D array depending on just 
%latitude and longitude 
%
%see also: getlev_ccm

function [x] = getlev(aray,lev);

%Declare global variables to be used in this function; must run script
%'rdccm_short' first
global XAX YAX l_pres ps;

aray=squeeze(aray);

%These are coefficients CCM3 uses to convert from sigma levels to 
%pressure levels
psa = [4.8093 13.0731 32.5591 63.9471 81.6768 78.0201 73.3671 67.6476 60.8624 53.1095 44.5995 35.6607 26.7266 18.3069 10.9421 5.1470 1.3519 0];
psb = [0 0 0 0 17.3664 60.6928 115.8237 183.5918 263.9851 355.8459 456.6760 562.5875 668.4428 768.2030 855.4653 924.1285 969.0938 992.5282];

%Arrays to hold the value of the variable at the levels immediately above
%and below the selected pressure level
x_up = zeros(length(YAX),length(XAX));
x_lo = zeros(length(YAX),length(XAX));


%Arrays to hold the pressure at the midpoint of the levels immediately
%above and below the selected pressure level
p_up = zeros(length(YAX),length(XAX));
p_lo = zeros(length(YAX),length(XAX));

%If selected pressure level is within the top 4 sigma levels, which depend
%only on pressure, x_up, x_lo, p_up, and p_lo can be assigned for all
%gridpoints at once
if (lev <= psa(4))
  if (lev > psa(3))
    x_up = squeeze(aray(3,:,:));
    x_lo = squeeze(aray(4,:,:));
    p_up = squeeze(l_pres(3,:,:));
    p_lo = squeeze(l_pres(4,:,:));
  elseif (lev > psa(2))
    x_up = squeeze(aray(2,:,:));
    x_lo = squeeze(aray(3,:,:));
    p_up = squeeze(l_pres(2,:,:));
    p_lo = squeeze(l_pres(3,:,:));
  elseif (lev >= psa(1))
    x_up = squeeze(aray(1,:,:));
    x_lo = squeeze(aray(2,:,:));
    p_up = squeeze(l_pres(1,:,:));
    p_lo = squeeze(l_pres(2,:,:));
  end

%Otherwise, we have to figure out which sigma levels the selected pressure
%level falls between for each gridpoint
else

%psx is the value that the surface pressure would have to be in order for
%the selected pressure level to exactly coincide with a given sigma level
  psx = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
  psx(5:18) = 1e5*(lev - psa(5:18))./psb(5:18);

%xps is 1 at a gridpoint if the selected pressure level falls between 
%the current pair of sigma levels, and 0 if it does not; pairs of sigma
%levels, beginning with 4 and 5 and ending witih 17 and 18,  are looped 
%through, with x_up, x_lo, p_up, and p_lo getting assigned values only 
%when the pressure level falls between the current pair of sigma levels
  xps = ps >= psx(5);
  x_up = x_up + xps.*squeeze(aray(4,:,:));
  x_lo = x_lo + xps.*squeeze(aray(5,:,:));
  p_up = p_up + xps.*squeeze(l_pres(4,:,:));
  p_lo = p_lo + xps.*squeeze(l_pres(5,:,:));
  for l = 6:18
    xps = (ps >= psx(l)) & (ps < psx(l-1));
    x_up = x_up + xps.*squeeze(aray(l-1,:,:));
    x_lo = x_lo + xps.*squeeze(aray(l,:,:));
    p_up = p_up + xps.*squeeze(l_pres(l-1,:,:));
    p_lo = p_lo + xps.*squeeze(l_pres(l,:,:));
  end
end

%Calculate the value of the variable interpolated to the selected pressure
%level for every gridpoint all at once, using the values of the variable at
%the sigma levels above and below the pressure level as well as the pressures
%at the sigma levels above and below the pressure level for each gridpoint
x = x_up + (log(lev*100) - log(p_up)).*(x_lo - x_up)./(log(p_lo) - log(p_up));
