function [xn,yn,zn]=advectp(X,Y,Z,z_coord,U_,V_,w,dt,x,y,z)

%keyboard

%xy2ll=0.01/1426.1; %conversion factor from lat lon to x y for middle of Monterey Bay
xy2ll=0.01/(m_lldist(x+[0;0.01],y+[0;0.01])*1e3);

sizeZ=size(Z);

%z_coord=nan(sizeZ(1),1);
%for zi=1:sizeZ(1)
%  z_coord(zi,1)=interp2(X,Y,squeeze(Z(zi,:,:)),x,y);
%end

idxnnan=~isnan(z_coord);
nl=interp1(z_coord(idxnnan),[1:sum(idxnnan)],z);%level number
if isnan(nl);
  nl=interp1(z_coord(idxnnan),[1:sum(idxnnan)],z,'nearest','extrap');
end
h=nl-floor(nl);%vertical interp coef

if nl==sizeZ(end)|nl==1
  u=interp2(X,Y,squeeze(U_(:,:,nl))',x,y);

  v=interp2(X,Y,squeeze(V_(:,:,nl))',x,y);

else
  u1=interp2(X,Y,squeeze(U_(:,:,floor(nl)))',x,y);

  u2=interp2(X,Y,squeeze(U_(:,:,ceil(nl)))',x,y);

  u=(1-h)*u1+h*u2;

  v1=interp2(X,Y,squeeze(V_(:,:,floor(nl)))',x,y);

  v2=interp2(X,Y,squeeze(V_(:,:,ceil(nl)))',x,y);

  v=(1-h)*v1+h*v2;
end

xn=u*dt*xy2ll+x;
yn=v*dt*xy2ll+y;
zn=w*dt+z;


