function FillIndex = roms_fillindex(data)
% FillIndex = roms_fillindex(DATA)
%
% For DATA on a 2-dimensional grid there may be missing values (NaN) 
% because of missing data in the original input. These need to be filled 
% before interpolation to the ROMS grid or there is potential for NaNs
% entering the ROMS values near the land/sea mask. This function will find 
% the 1-D indicies such that DATA = DATA(FillIndex) will fill the entries 
% in DATA that are NaN using a griddata 'nearest' method applied to the 
% good data.
%
% The FillIndex values can be computed once and saved for speedy filling.
%
% To see how it works, just pcolor the FillIndex result and you'll get the
% idea
%
% John Wilkin 
% Apr 2008

if ndims(data)~=2
  error('Input must be 2-dimensional')
end
I = 1:numel(data);
I = reshape(I,size(data));
good = find(isnan(data)~=1);
bad = find(isnan(data)==1);
[nx ny] = size(data);
[X Y] = meshgrid(1:nx,1:ny);
X=X'; Y=Y';
I(bad) = griddata(X(good),Y(good),I(good),X(bad),Y(bad),...
  'nearest',{'QJ'});
FillIndex = I;
