function [Tm,xc,yc,zc] = interp3ddata(xl,yl,zl,TL,xc,yc,zc,msk);
%
% [T,x,y,z] = interp3ddata(xl,yl,zl,TL,xc,yc,zc,msk);
% dataname is a string ('temp' or 'salt')

[nxl,nyl,nzl,nt]=size(TL);
nxc=size(xc,1);
nyc=size(yc,2);
nzc=prod(size(zc));

% Use NaN's to represent no data
TL=sq(TL);

% $$$ % interpolate vertically first
% $$$ disp('Vertical interpolation')
% $$$ TLL=reshape(TL,[nxl*nyl nzl]);
% $$$ TN = zeros(nxl*nyl,nzc);
% $$$ TL = TN*0;
% $$$ for i=1:size(TLL,1)
% $$$   iz = find(~isnan(TLL(i,:)));
% $$$   if length(iz)>1
% $$$     TL(i,:)=interp1(zl(iz),TLL(i,iz)',zc,'linear');
% $$$     TN(i,:)=interp1(zl(iz),TLL(i,iz)',zc,'nearest','extrap');
% $$$   end
% $$$ end
% $$$ in = find(isnan(TL)); TL(in) = TN(in);
% $$$ TL=sq(reshape(TL,[nxl nyl nzc]));
% $$$ clear TLL TN



% Create halo in X direction
disp('Adding wrap around data')
xl=[xl(end-2:end)-360 xl xl(1:3)+360];
for k=1:size(TL,3);
 TLL(:,:,k)=[TL(end-2:end,:,k)' TL(:,:,k)' TL(1:3,:,k)']';
end
TL=TLL; clear TLL


% Interpolate horiztonally
%load MASKS

for k=1:size(TL,3);
 disp( sprintf('Interpolating (x-y) level %i',k) );
 Tl=interp2(yl,xl,TL(:,:,k),yc,xc,'cubic');
 Te=interp2(yl,xl,TL(:,:,k),yc,xc,'linear');
 jj=find( isnan(Tl) ) ; Tl(jj)=Te(jj);
 Te=interp2(yl,xl,TL(:,:,k),yc,xc,'nearest');
 jj=find( isnan(Tl) ) ; Tl(jj)=Te(jj);
 %Tl(find(isnan(Tl)))=0;
 Tl=xyexpand(Tl.*sq(msk(:,:,1)),11);
 Tm(:,:,k)=Tl.*sq(msk(:,:,1));
end
clear Tl Te jj TL

%load VGRID.mat nzc zc 

disp('Vertical interpolation')
Tl=reshape(Tm,[nxc*nyc nzl]);
Tm=interp1(zl,Tl',zc,'linear','extrap');
Tn=NaN*Tm;
for i=1:size(Tl,1)
  iz = find(~isnan(Tl(i,:)));
  if length(iz)> 1
    Tn(:,i)=interp1(zl(iz),Tl(i,iz)',zc,'nearest','extrap');
  end
end
in = find(isnan(Tm)); Tm(in) = Tn(in);
Tm=reshape(Tm',[nxc nyc nzc]);
Tm=Tm.*msk;
clear Tl Tn

Tmsk=mask(Tm);
it=0;
while( sum(Tmsk(:)-msk(:)) & it<20)
it=it+1;
  disp( sprintf('Number of points needing vertical extrapolation = %i', ...
		sum(msk(:)-Tmsk(:)) ) );
  Te=2*Tm(:,:,[1 1:end-1])-Tm(:,:,[1 1 1:end-2]);
  jj=find( Tmsk==0 & msk==1 );
  Tm(jj)=Te(jj);
  Tmsk=mask(Tm);
end

  
% Remove NaN's
Tm( find(isnan(Tm)) )=0;