function [Tm,xc,yc,zc] = interp3ddata_int(xl,yl,zl,TL,xc,yc,zc,msk);
%
% [T,x,y,z] = interp3ddata_int(xl,yl,zl,TL,xc,yc,zc,msk);

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

% Back out the grid spacing for model
dz(1)=-zc(1)*2;
zf(1)=0;
for k=1:nzc-1;
 zf(k+1)=zf(k)-dz(k);
 dz(k+1)=2*(zf(k+1)-zc(k+1));
end
zf(nzc+1)=zf(nzc)-dz(nzc);

% Assume some grid spacing for Levitus
zlf=[0 (zl+zl([2:end end]))/2];
dzl=zlf(1:nzl)-zlf(2:nzl+1);

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

% 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
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


disp('Vertical interpolation')
Tm( find(isnan(Tm)) )=0;
hm=reshape( -sum(reshape(msk,[nxc*nyc nzc])*dz',3) ,[nxc nyc]);
hl=reshape( -sum(reshape(mask(Tm),[nxc*nyc nzl])*dzl',3) ,[nxc nyc]);
for kk=1:nzc;
 z1m=max(zf(kk),hm);
 z2m=max(zf(kk+1),hm);
 dm=max(z1m-z2m,0);
 dw=0;
 ti=0;
 for k=1:nzl;
  z1l=max(zlf(k),hl);
  z2l=max(zlf(k+1),hl);
  z1=min(z1m,z1l);
%  z2=min(z2m,z2l);
  z2=max(z2m,z2l);
  dd=max(z1-z2,0);
  dw=dw+dd;
  ti=ti+dd.*Tm(:,:,k);
 end
 Tl(:,:,kk)=ti./sq(dw);
end
Tm=Tl.*msk;
clear Tl

Tmsk=mask(Tm);
it=0;
while( sum(Tmsk(:)-msk(:)) & it<15)
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;