% out = convert_sigma_to_pres(in, ps, newlevs);
% 
%   This function converts CCM data from hybrid-sigma coordinates 
%   to sigma coordinates.  Parameters hyam, hybm, PO are
%   predefined in the script.
%   from Daniel J. Vimont
%   http://www.aos.wisc.edu/~dvimont/matlab/Utility_Tools/convert_sigma_to_pres.html
% 
%   made changes by volley to suit for SCAM output for dimensions.
%   2014.10.7
% 
%   Input Variables:
%   in = data to be converted
%   ps = surface pressure data (CAREFUL:  we will attempt to convert this
%        to hPa from Pa).
%   newlevs = levels to interpolate to, in hPa
% 
%   Output Variables:
%   out = data on sigma levels.
% 
%   Conventions are:
%   1.  Assume ps is in Pa (we will convert this to hPa)
%   2.  Assume newlev is in hPa
%   3.  [nlon, nlat, nlev, ntim] = size(in);
%   4.  [nlon, nlat, ntim] = size(ps);
%   5.  nlev2 = length(newlevs);
%   6.  [nlon, nlat, nlev, ntim] = size(out);
%   
%   Note that for the above, singleton dimensions should be preserved.
% 
% Listing of function convert_sigma_to_pres

function out = convert_sigma_to_pres(in, ps, newlevs, hyam, hybm)

% hyam = [0.00480929994955659, 0.013073099777102, 0.032559100538492, ...
% 	0.0639470964670177, 0.0802583247423168, 0.0766648948192596, ...
% 	0.0720923840999599, 0.0664718076586719, 0.0598041415214539, ...
% 	0.0521853193640709, 0.0438226312398908, 0.035038441419601, ...
% 	0.0262589752674099, 0.017985042184591, 0.010747664608061, ...
% 	0.00505276303738358, 0.0013234377838671, 0]';
% 
% 
% hybm = [0, 0, 0, 0, 0.018784884363413, 0.062048017978668, ...
% 	0.11709871888161, 0.18476781249046, 0.265043288469309, ...
% 	0.356770157814029, 0.457452863454819, 0.56321024894714, ...
% 	0.668910741806027, 0.768524885177609, 0.855659365653988, ...
% 	0.92422324419022, 0.96912252902985, 0.992528021335598]';

pnot = 1000;

nlev = length(hyam);
newlevs = newlevs(:);
nlev2 = length(newlevs);

%  Go through argument check.
size_ps = size(ps);
size_in = size(in);

if ndims(ps) ~= (ndims(in)-1) 
  error('Variable ''in'' should have one more dimension than ''ps''.')
end

if size_in(3) ~= nlev;
  error(['Either ''in'' has an inconsistent level dimension, or the' ...
	 'order (nlon, nlat, nlev, ntim) is wrong.                 '])
end

if length(size_in)==3
    size_in(4)=1;
%     disp(['the time dimension might be 1, if it is not, please check']);
end

%  Reshape so level is the first dimension
in = shiftdim(in, 2); 
in = reshape(in, size_in(3), prod(size_in)/nlev);

ps = shiftdim(ps, 2);
ps = reshape(ps, 1, prod(size_ps));

%  Check for size consistency
if size(ps, 2) ~= size(in, 2);
  error('Variables ''in'' and ''ps'' have inconsistent dimensions')
end

npts = prod(size_in)/nlev;

%  Now that ps is 2D, make sure units are in hPa:
if round(log10(mean(mean(ps)))) == 5;
  ps = ps/100;
end

%  Now, do interpolation
out = repmat(NaN, [nlev2 prod(size_in)/size_in(3)]);
plevs = pnot*hyam*ones(1, npts)+hybm*ps;

%  This method is considerably faster than the following - both produce
%  the same results.
for i = 1:nlev2;
  
  for j = 1:npts
    
    xup = sum(plevs(:,j) < newlevs(i));
    
    if xup < nlev;
      out(i,j) = in(xup,j) + ...
	  (log(newlevs(i)) - log(plevs(xup,j))) * ...
	  (in(xup+1,j) - in(xup,j))/(log(plevs(xup+1,j))-log(plevs(xup,j)));
    end
    
  end
  
end

%  This method is considerably slower than the preceding
% for i = 1:npts
%  out2(:, i) = interp1(log(plevs(:,i)), in(:,i), log(newlevs));
% end

%  Reshape back
out = reshape(out, [nlev2 size_in(4) size_in(1:2)]);
% for nlat and nlon are not 1
% out = shiftdim(out, dims_in - 2);
% for nlat=nlon=1
out = permute(out,[3 4 1 2]);
