function [i,j,z]=rcloseto(fn,tlat,tlon)
% RCLOSETO - Finds index of lat/lon in grd.nc file
% [i,j,z]=rcloseto(fn,tlat,tlon)
nc = netcdf(fn);
lat = nc{'lat_rho'}(:,:);
lon = nc{'lon_rho'}(:,:);
h =  nc{'h'}(:,:);

% calculate distance using algorithm in earthdist.m
radius = 6371*1000;
RCF = 180 / pi;
tlonr = tlon / RCF;
tlatr = tlat / RCF;
lonr = lon / RCF;
latr = lat / RCF;

c = cos(tlatr);
ax = c .* cos(tlonr);
ay = c .* sin(tlonr);
az = sin(tlatr);

c = cos(latr);
bx = c .* cos(lonr);
by = c .* sin(lonr);
bz = sin(latr);

ddist = acos(ax.*bx + ay.*by + az.*bz) .* radius;

[i,j]=find( ddist == min(ddist(:) ) );
for ii = 1:length(i),
  for jj = 1:length(j),
fprintf(1,'i = %d, target: %f ROMS: %f\n',i(ii),tlat,lat(i(ii),j(jj)));
fprintf(1,'j = %d, target: %f ROMS: %f\n',j(jj),tlon,lon(i(ii),j(jj)));
fprintf(1,'Depth of ROMS grid cell: %f\n',h(i(ii),j(jj)));
  end
end
z = h(i,j);
close(nc)