function dataout = roms_zinterp_to_newbathy(datain,oldh,grd_new,scoord)
% Interpolate a roms variable from one bathymetry to another, or change the
% scoordinate stretching parameters
%
% Assumes the horizontal grid doesn't alter
% just need to (a) do the vertical interpolation to new s-levels
%              (b) handle possible extrapolation where hnew>hold
%
% Load inputs values into 'datain' and 'zin'
% interpolate to the new depths 'z' and store in 'data'
% [Don't need to interpolate the 2-D variables zeta, ubar, vbar]
%
% John Wilkin


% Create a new grid structure that computes z_r etc using the old
% bathymetry and/or new s-coordinate stretching parameters (i.e. this
% function can be used to convert data on the same grid and depth but
% to a new set of stretching values or number of vertical levels
grd_old = grd_new;

if ~isempty(oldh)
  % Take 'old' bathymetry from the inputs
  grd_old.h = oldh;
end

if nargin < 4
  % use the vertical stretching in the input grd structure
  scoord = [grd_new.theta_s grd_new.theta_b grd_new.Tcline grd_new.N];
end

% compute z_r etc
grd_old = roms_get_grid(grd_old,scoord,0,1);

% determine c-grid position of the input datain
pos = roms_cgridpos(datain,grd_old);
m   = getfield(grd_old,['mask_' pos]);
zin = getfield(grd_old,['z_' pos(1)]);
zin = reshape(zin,[size(zin,1) length(zin(1,:))]);
zout = getfield(grd_new,['z_' pos(1)]);
shapeout = size(zout);
zout = reshape(zout,[size(zout,1) length(zout(1,:))]);

% copy top and bottom rows of input data and input z so that if there are
% new z values out of range of the inputs, the interpolation does not 
% return NaN. What it will give is an upward and downward copy of the 
% shallowest and deepest available data (to limits -10000 m and +10 m
zin = [-10000*ones([1 length(zin)]); zin; 10*ones([1 length(zin)])];
datain = datain([1 1:end end],:);

% for speed, only compute in the water
water = find(m~=0);
dataout = zeros(size(zout));
for i=water'
  dataout(:,i) = interp1(zin(:,i),datain(:,i),zout(:,i));
  if rem(i,1000)==0
    disp([ ' On ' int2str(i) ' of ' int2str(length(datain)) ...
      ' points changing z stretching for new bathy' ])
  end
end
dataout = reshape(dataout,shapeout);

% check for NaNs
bad = find(isnan(dataout)==1);
if ~isempty(bad)
  warning([' There are ' int2str(length(bad)) ' NaNs in the output'])
end
