function VARinterp=interpTPXO(VAR,harmonic,TPXOfileVAR,TPXOfileGRID,lon,lat,mask)
%*****************************************************************
% Extract TPXO data and interpolate onto ROMS grid
% where possible. TriScatteredInterp will produce zeros or NaNs
% if unable to interpolate (it will not extrapolate up into bays etc.) -
% in those place set the values equal to the nearest nonzero value.
% Inputs:
% VAR: hRE, hIm, uRE, uIm, vRe, or vIm
% harmonic: one of: M2 S2 N2 K2 K1 O1 P1 Q1 MF MM M4 MS4 MN4
% lon, lat: 2D arrays of longitude and latitude to interpolate to
% mask: array of 1s and 0s corresponding to wet and dry points of lon & lat
% Output: VARinterp (VAR on lon, lat grid)
% e.g. zr=extract_tpxo('hRe','M2','h_tpxo72.nc','grid_tpxo71.nc',lon,lat,rmask); %interp Re part of M2 elevation
%        ur=extract_tpxo('uIm','S2','u_tpxo72.nc','grid_tpxo71.nc',lon,lat,rmask);  %interp Im part of S2 u-velocity
% Untested for domains extending across dateline - probably bogus there.
% TPXO files must be on the matlab path
% J. Luick November 2011
%*****************************************************************

TPXOnames={'M2' 'S2' 'N2' 'K2' 'K1' 'O1' 'P1' 'Q1' 'MF' 'MM' 'M4' 'MS4' 'MN4'};
Nharmonic=find(strcmp(harmonic,TPXOnames));

lonmin=min(min(lon));
lonmax=max(max(lon));
latmin=min(min(lat));
latmax=max(max(lat));

VARinterp=mask;

%Extract TPXO data
%Step 1 - interpolate where possible
if strcmp(VAR,'hRe') ||  strcmp(VAR,'hIm')     %real or imag part of complex elevation "h"
   X=nc_varget(TPXOfileVAR,'lon_z');
   Y=nc_varget(TPXOfileVAR,'lat_z');
   Z=nc_varget(TPXOfileVAR,VAR);
   Z=squeeze(Z(Nharmonic,:,:));
   ff=find(Y>=latmin-0.5 & Y<=latmax+0.5 & X>=lonmin-0.5 & X<=lonmax+0.5);
   X=X(ff);
   Y=Y(ff);
   Z=Z(ff);
   F=TriScatteredInterp(X,Y,Z);    
   VARinterp=F(lon,lat);    %z=F(x,y) interpolates to 1-d vectors x and y  (z also is 1-d) 
elseif strcmp(VAR,'URe') ||  strcmp(VAR,'UIm') %real part of complex u-component of 2D velocity
   X=nc_varget(TPXOfileVAR,'lon_u');
   Y=nc_varget(TPXOfileVAR,'lat_u');
   Z=nc_varget(TPXOfileVAR,VAR);
   H=nc_varget(TPXOfileGRID,'hu');
   Z=squeeze(Z(Nharmonic,:,:));
   ff=find(Y>=latmin-0.5 & Y<=latmax+0.5 & X>=lonmin-0.5 & X<=lonmax+0.5 & H>0);
   X=X(ff);
   Y=Y(ff);
   Z=Z(ff)./H(ff); %Convert transport to 2D velocity
   F=TriScatteredInterp(X,Y,Z);    
   VARinterp=F(lon,lat);    
elseif strcmp(VAR,'VRe') ||  strcmp(VAR,'VIm') %real part of complex u-component of 2D velocity
   X=nc_varget(TPXOfileVAR,'lon_v');
   Y=nc_varget(TPXOfileVAR,'lat_v');
   Z=nc_varget(TPXOfileVAR,VAR);
   H=nc_varget(TPXOfileGRID,'hv');
   Z=squeeze(Z(Nharmonic,:,:));
   ff=find(Y>=latmin-0.5 & Y<=latmax+0.5 & X>=lonmin-0.5 & X<=lonmax+0.5 & H>0);
   X=X(ff);
   Y=Y(ff);
   Z=Z(ff)./H(ff);
   F=TriScatteredInterp(X,Y,Z);    
   VARinterp=F(lon,lat);    
end
%Step 2 - elsewhere use nearest good value
ff=find (mask==1 & VARinterp==0);  %ff: wet points that were not assigned 0
for n=1:length(ff);
   dr=spheric_dist(lat(ff(n)),lat,lon(ff(n)),lon);
   [j,i]=find(dr==(min(min(dr(VARinterp~=0)))));
   VARinterp(ff(n))=VARinterp(j,i);
end
ff=find (mask==1 & isnan(VARinterp));  %ff: wet points that were assigned NaN
for n=1:length(ff);
   dr=spheric_dist(lat(ff(n)),lat,lon(ff(n)),lon);
   [j,i]=find(dr==(min(min(dr(~isnan(VARinterp))))));
   VARinterp(ff(n))=VARinterp(j,i);
end
VARinterp(isnan(VARinterp))=0;

return