function vgeo_grd(yb, fn, x, grd, nf, plt, inst, info)
%------------------------------------------------------------------
% NAME
%
%   vgeo_grd - Plot data from a grd vgeo file
%
% SYNOPSIS
%
%   vgeo_grd(year_base, file_name, x, grd, filter, title, institute)
%
% PARAMETERS
%
%   year_base
%
%   file_name
%      Name of vgeo grd file, i.e. output of programme vgeo when grids are
%      specified. This file has the following columns:
%
%         1. time (decimal day)
%         2. longitude
%         3. latitude
%         4. depth
%         5. geostrophic current
%
%      Only give the root name, since suffix '.grd' is appended inside the
%      function. 
%
%   x
%      Index to x variable, according to the columns in the data file (i.e.
%      most probably 1, 2 or 3). If negative, the axis will be reversed.
%     
%   grd
%      Vector with grid values to plot. If empty, a list with all different
%      depths in the data file is generated.
%
%   filter
%      A number (>=2) for filtering data.
%
%   title
%
%   institute
%
%   info
%      Information string.
%
% DESCRIPTION
%
%   Generate a plot of geostrophic current at the specified depths (all of them
%   if no depths are specified). Velocities are first averaged together to
%   remove large variations, and lowpass filtered to remove high frequency
%   oscillations. 
%
% OUTPUT
%
%   <file_name>_grd.ps
%      Output postscript file.
%
%------------------------------------------------------------------

lines   = { '-', '-.', '--', ':' };
l_names = { 'solid', 'dashdot', 'dashed', 'dotted' };

reverse = 0;
if x < 0,
  reverse = 1;
  x = abs(x);
end
ps_file = [ fn '_grd.ps' ];

%--- Read data:
file = [ fn '.grd' ];
d = read_asc(file);
d = tonan(d);
nc = size(d, 2);
if nc ~= 5,
	error([ 'Wrong number of columns in file ' file ]);
end

%--- Indexes to array:
t = 1; z = 4; v = 5;
d = d';

%--- Select data:
if ~length(grd),
	%--- Get the different depth values:
	i = 1;
	grd = -1;
	zz = d(z,i);
	while ~length(find(grd == zz)),
		grd = [ grd, zz ];
		i   = i+1;
		zz  = d(z,i);
	end
  grd = grd(2:length(grd));	
end

%--- Create figure:
xdir_str = 'normal';
if reverse,
  xdir_str = 'reverse';
end
figure;
orient landscape;
axes('position', [0.1, 0.1, 0.8, 0.7], 'XDir', xdir_str);
hold on;

%--- Plot velocity:
lgstr = '';
ii = 1:length(grd);
for i=ii,
  j = mod(i,length(lines)+1);
  if j == 0,
    j = 1;
  end
	ls = lines{j};
	iz = find(d(z,:) == grd(i));
	%--- Calculate running means:
	[vv, good] = rmean(d(v,iz), nf);
	xx = rmean(d(x,iz(good)), nf);
	%--- Filter data:
  if nf > 1,
    [b,a] = ellip(4, 0.1, 50, 2/nf);
    vf = filter(b,a, vv);
  else
    vf = vv;
  end
	%--- Plot data:
	plot(xx, vf, ls)
	%--- Upadate legend string:
	lgstr = [ lgstr ' - [' l_names{j} ': ' sprintf('%-5.2d',grd(i)) 'm]' ];
end
switch(x)
	case 1,
		xlabel([ 'Decimal Days (' int2str(yb) ')' ]);
	case 2,
		xlabel('Longitude');
	case 3,
		xlabel('Latitude');
end
grid on

%--- Get name of output file:
%ps_file = next_psf([ fn '.ps' ]);
lgstr = [ lgstr '  - PS File: ' ps_file ];
info = [ 'Geostrophic Current - ' info ];

%--- Plot residual things:
subtitle = [ to_date(yb, d(t,1)) '  to  ' ...
			to_date(yb, d(t,length(d(t,:)))) ];
plotstuf('-title', plt, '-subtitle', subtitle, '-note', lgstr, ...
		'-institute', inst, '-info', info, '-format', 'landscape');

%--- Print output:
eval([ 'print -dps ' ps_file ]);




