function [val]=idma_interp_2d(lon,lat,fld,method);
%IDMA_INTERP_2D linearly interpolates gcmfaces field (fld) to
%   a set of locations (lon, lat) using one of several methods:
%   'natural', 'linear', 'nearest', or 'mix' (default).
%   The 'mix' is `natural' extended with 'nearest' when
%   the input field has been land-masked with NaNs.
%
%Example:
%     lon=[-179.9:0.2:179.9]; lat=[-89.9:0.2:89.9];
%     [lat,lon] = meshgrid(lat,lon);
%     fld=mygrid.Depth.*mygrid.mskC(:,:,1);
%     [val]=idma_interp_2d(lon,lat,fld);
%     figureL; pcolor(lon,lat,val); shading flat;

gcmfaces_global;
if isempty(whos('method'));
    method='mix';
end;

if isempty(which('DelaunayTri'));
    error('this code needs DelaunayTri that is not found');
end;

XC=convert2array(mygrid.XC);
YC=convert2array(mygrid.YC);
VEC=convert2array(fld);
val=NaN*lon;
for ii=1:3;
    if ii==1;
        myXC=XC; myXC(myXC<0)=myXC(myXC<0)+360;
        jj=find(lon>90);
    elseif ii==2;
        myXC=XC;
        jj=find(lon>=-90&lon<=90);
    else;
        myXC=XC; myXC(myXC>0)=myXC(myXC>0)-360;
        jj=find(lon<-90);
    end;
    kk=find(~isnan(myXC));
    TRI=DelaunayTri(myXC(kk),YC(kk));
    if strcmp(method,'mix');
        F = TriScatteredInterp(TRI, VEC(kk),'natural');
        tmp1=F(lon(jj),lat(jj));
        F = TriScatteredInterp(TRI, VEC(kk),'nearest');
        tmp2=F(lon(jj),lat(jj));
        tmp1(isnan(tmp1))=tmp2(isnan(tmp1));
        val(jj)=tmp1;
    else;
        F = TriScatteredInterp(TRI, VEC(kk),method);
        val(jj)=F(lon(jj),lat(jj));
    end;
end;