function sgeo_grd(vg_file, out_file, yb, x, s0, v_scale, s_scale, ...
    v_dir, depths, plt, info, inst)
%  1.  NAME
%
%  sgeo_grd - Calculate and plot geostrophic slope height
%
%  2.  SYNOPSIS
%
%  sgeo_grd(vg_file, out_file, yb, x, s0, v_scale , s_scale, v_dir,
%  depths, plt, info, inst )
%
%  3.  DESCRIPTION
%
%  Read geostrophic current data from a grid ascii file, as output by C
%  utility vgeo(1) and compute the corresponding slope heights. Generate
%  a plot for each of the specified depths and save data into an ascii
%  file.
%
%  The input file must contain the following columns:
%
%  1. Time in decimal days
%
%  2. Longitude in degrees
%
%  3. Latitide in degrees
%
%  4. Depth in meters
%
%  5. Geostrophic current
%
%  This function requires a lot of parameters. The best way is to use it
%  from Matlab script sgeo(1m).
%
%  4.  PARAMETERS
%
%     vg_file
%        Name if input file with geostrophic current data.
%
%     out_file
%        Root name of output file. Suffixes are automatically added.
%
%     yb Base year used to determine the decimal days in the input file.
%
%     x  Column number to x-axis for plotting. It must be one of 1,2 or
%        3. If negative,the axis direction will be reversed.
%
%     s0 Initial slope value in meters.
%
%     v_scale
%        Factor to get geostrophic currents from input file in m/s.
%
%     s_scale
%        Scaling factor for plotting slope height points.
%
%     v_dir
%        The direction that was used to compute the geostrophic current.
%        This is necessary in order to get the correct slope variations.
%        This is one of northward,southward,westward of eastward.
%     depths
%        Array with depths values at which plotting of slope heights is
%        desired.  The values must be present in the input file (i.e. no
%        interpolation is carried out here). If empty,all available
%        depths in the input file will be plotted.
%
%     plt
%        Title of plot.
%
%     info
%        Additionnal information string.
%
%     inst
%        Institute stamp.
%
%  5.  OUTPUT
%
%  The following files are generated:
%
%     out_file.ps
%        Postcript output of plot.
%
%     out_file.asc
%        Ascii data of slope height values.
%
%  6.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

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

reverse = 0;
if x < 0,
  reverse = 1;
  x = abs(x);
end
ps_file = [ out_file '.ps'  ];
asc_file= [ out_file '.asc' ];

%--- Read the data:
d = tonan(read_asc(vg_file));
d = d';
nl = size(d,2);

%--- Indexes:
t = 1; lon = 2; lat = 3; z = 4; v = 5;

%--- Scale velocities:
d(v,:) = d(v,:) * v_scale;

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

%--- Prepare ascii output:
fid = fopen(asc_file, 'w');
fprintf(fid, '%s\n', ...
    [ '%--- sgeo.m: Calculate slope from geostrophic current.' ], ...
    [ '%    1. Time in decimal days (', int2str(yb), ')' ], ...
    [ '%    2. Longitude' ], ...
    [ '%    3. Latitude' ], ...
    [ '%    4. Depth' ], ...
    [ '%    5. Geostrophic Current (m/s)' ], ...
    [ '%    6. Distance (m)' ], ...
    [ '%    7. Slope Height (m)' ], ...
    [ '%--- sgeo.m: ', datestr(now) ]);
fprintf(fid, '\n');

%--- Create figure:
figure;
orient tall;

orient landscape;
xdir_str = 'normal';
if reverse,
  xdir_str = 'reverse';
end
axes('position', [0.1, 0.1, 0.8, 0.7], 'XDir', xdir_str);
hold on;
grid on;

%--- Loop through the specified depths:
lgstr = '';
ii = 1:length(depths);
for i=ii,
  
  %--- Select correct data:
  iz = find(d(z,:) == depths(i));
  if isempty(iz),
    error([ 'No data for depth ' sprintf('%-f', depths(i)) ]);
  end
  
  %--- Sort indexes to get increasing x-values:
  [xs, ix] = sort(d(x,iz));
  iz = iz(ix);
  
  %--- Select only good values:
  jj = find(~isnan(d(v,iz)));
  iz = iz(jj);
  
  %--- Eventually, reverse the indexes so that integration is carried out 
  %    in the same direction as the derivative of the geostrophic current.
  %
  % NOTE: ds has to change sign depending on the direction. This is because
  %       geostrophic caluclations with vgeo produce positive values for
  %       currents in the direction startboard to the specified track. The
  %       change of sign is necessary to get the correct slope changes.
  switch(v_dir)
    case 'southward',
      if d(lat,1) < d(lat,nl),
        iz = iz( length(iz):-1:1 );
      end
      ds = -1;
    case 'northward',
      if d(lat,1) > d(lat,nl),
        iz = iz( length(iz):-1:1 );
      end
      ds = 1;
    case 'eastward',
      if d(lon,1) > d(lon,nl),
        iz = iz( length(iz):-1:1 );
      end
      ds = -1;
    case 'westward',
      if d(lon,1) < d(lon,nl),
        iz = iz( length(iz):-1:1 );
      end
      ds = 1;
  end
  
  %--- Calculate distance between each points:
  nz = length(iz);
  ii = 2:nz;
  alat = 0.5 * (d(lat,iz(ii)) + d(lat,iz(ii-1)));
  dlon = d(lon,iz(ii)) - d(lon,iz(ii-1));
  dlat = d(lat,iz(ii)) - d(lat,iz(ii-1));
  dx = lon_to_m(dlon, alat);
  dy = lat_to_m(dlat, alat);
  ds = ds * sqrt(dx .* dx + dy .* dy);
  ds = [ 0 ds ];

  %--- Integrate the data:
  s = ones(1,nz) * NaN;
  s(1) = 0;
  for j=ii,
    s(j) = s(j-1) + 0.5 * (d(v,iz(j)) + d(v,iz(j-1))) * ds(j);
  end
  s = s .* coriolis(d(lat,iz)) ./ gravity(d(z,iz), d(lat,iz)); 
  
  %--- Set the vertical position of the slope:
  s = s + s0;
  
  %--- Scale slope data:
  s = s * s_scale;

  %--- Plot the slope:
  j = mod(i,length(lines)+1);
  if j == 0,
    j = 1;
  end
  ls = lines{j};
  plot(d(x,iz), s, ls);
  
  %--- Update legend string:
  lgstr = [ lgstr  ' - [' l_names{j} ': ' sprintf('%-5.2d', depths(i)) 'm]' ];
  
  %--- Save ascii data:
  jj=1:length(s);
  for j=jj,
    fprintf(fid,'%15.6f ', d(t,iz(j)), d(lon,iz(j)), d(lat,iz(j)), ...
        depths(i), d(v,iz(j)), ds(j), s(j)/s_scale);
    fprintf(fid, '\n');
  end
  
end
fclose(fid);

switch(x)
	case 1,
		xlabel([ 'Decimal Days (' int2str(yb) ')' ]);
	case 2,
		xlabel('Longitude');
	case 3,
		xlabel('Latitude');
end

%--- Get name of output file:
lgstr = [ lgstr '  - PS File: ' ps_file ];
if length(info) > 0,
  info = [ 'Geostrophic Slope - ' info ];
end

%--- 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 ]);













