function seagrid2ecom(theSeagridFile, theEcomFile)

% seagrid2ecom -- Convert "seagrid" output to "ecom" input".
%  seagrid2ecom('theSeagridFile', 'theEcomFile') converts
%   the given SeaGrid file to an "ecom" grid file.  Where
%   filenames are not provided, the "uigetfile/uiputfile"
%   dialogs are invoked.  The data represent a grid with
%   twice as many grid lines as needed for ECOM, so that
%   we can work with cell-centers.
 
% Copyright (C) 1999 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 07-Jun-1999 15:43:58.
% Updated    26-Oct-1999 14:54:49.
% Updated    23-Mar-2000 11:12:15.  by RPS (rsignell@usgs.gov) 
%   to add appropriate "model_grid" header stuff, fix bogus mask subsampling,
%   fix reversed IM and JM, and handle remaining NaNs in the bathymetry.

RCF = 180 / pi;   % Radian conversion factor.

if nargin < 1, theSeagridFile = '*.mat'; end
if nargin < 2, theEcomFile = 'ecom_model_grid.txt'; end

% Get filenames.

if any(theSeagridFile == '*')
	[f, p] = uigetfile(theSeagridFile, 'Select SeaGrid File:');
	if ~any(f), return, end
	if p(end) ~= filesep, p(end+1) = filesep; end
	theSeagridFile = [p f]
end

if nargin < 2
	[f, p] = uiputfile(theEcomFile, 'Save to Ecom File:');
	if ~any(f), return, end
	if p(end) ~= filesep, p(end+1) = filesep; end
	theEcomFile = [p f]
end

disp([' ## SeaGrid Source File  : ' theSeagridFile])
disp([' ## Ecom Destination File: ' theEcomFile])

% Load input file.

try
	theSeagridData = load(theSeagridFile, 's');
catch
	disp([' ## Unable to load: "' theSeagridFile '"'])
	return
end

% Open output file.

fecom = fopen(theEcomFile, 'w');

if fecom < 0
	disp([' ## Unable to open: "' theEcomFile '"'])
	disp(' ## If already open elsewhere, please close, then try again.')
	return
end

% How rigid is the following format?
% For large grid cells, the H1 and H2 fields collide.

theFormat = ...
	'%4.0f %4.0f %10.2f %10.2f %10.2f %8.1f %8.1f %8.1f %15.6f %15.6f\n';

s = theSeagridData.s;

lon = s.geographic_grids{1};
lat = s.geographic_grids{2};

% Cell sizes.

dx = earthdist(lon(:, 1:2:end-1), lat(:, 1:2:end-1), lon(:, 3:2:end), lat(:, 3:2:end));
dy = earthdist(lon(1:2:end-1, :), lat(1:2:end-1, :), lon(3:2:end, :), lat(3:2:end, :));

h1 = dx;   % ECOM nomenclature.
h2 = dy;

% Cell orientations: need to check indices and directions here.

dlon = diff(lon.').';
dlat = diff(lat.').';
clat = cos(lat / RCF);
clat(:, end) = [];
ang = atan2(dlat, dlon .* clat) * RCF;
ang = 0.5 * (ang(2:2:end, 1:2:end) + ang(2:2:end, 2:2:end));

% Cell-centers.

lon = lon(2:2:end-1, 2:2:end-1);
lat = lat(2:2:end-1, 2:2:end-1);

% Bathymetry and mask.
bathymetry=s.gridded_bathymetry;

if isempty(bathymetry)
	bathymetry = zeros(size(lon));
end

mask = s.mask;
if isempty(mask)
	mask = zeros(size(bathymetry));
end

% Adjust sizes.
if all(size(bathymetry) > size(lon))
	bathymetry = bathymetry(2:2:end-1, 2:2:end-1);
end
if all(size(mask) > size(lon))
	mask = mask(2:2:end-1, 2:2:end-1);
end

min_depth = s.clipping_depths(1);
max_depth = s.clipping_depths(2);
bathymetry(bathymetry < min_depth) = min_depth;
bathymetry(bathymetry > max_depth) = max_depth;
bathymetry(logical(mask)) = -99999;

% RPS I don't know why there are still some nans in 
% the bathymetry at this point, since the mask should
% have caught them, but make'em land anyway.
bathymetry(isnan(bathymetry))=-99999;



% Coriolis parameter (use latitude).

coriolis = lat;

% Write the number of grid cells.

[rows, columns] = size(lon);

try
	fprintf(fecom, 'Created by SeaGrid on %s \n',date);
	fprintf(fecom, 'Vertical Segmentation - Sigma Levels (KB)\n');
	fprintf(fecom, '11\n');
	for sigma=-0:-0.1:-1;
          fprintf(fecom, '%4.2f\n',sigma);
        end
	fprintf(fecom, ...
            'Horizontal Segmentation (IM & JM) -- Grid Information\n');
	fprintf(fecom, '%5i%5i\n', columns, rows);
catch
	disp([' ## ' lasterr])
	fclose(fecom);
	return
end

% Be sure the -99999 stuff is correct.

% Output the grid.

for column = 1:columns   % Start at the left.
	remaining = column;
	if remaining == columns | (rem(remaining, 100) == 0)
		disp([' ## Columns remaining: ' int2str(remaining)])
	end
	for row = rows:-1:1   % Start at the bottom.
		try
			i = column;   % ECOM notation: i -> x; j -> y.
			j = rows-row+1;
			fprintf(fecom, theFormat, ...
	          	i, j, ...
			 	h1(row, column), h2(row, column), ...
			  	bathymetry(row, column), ...
			  	ang(row, column), ...
			  	coriolis(row, column), ...
			  	0.0, ...
			  	lon(row, column), lat(row, column));
		catch
			disp([' ## ' lasterr])
			fclose(fecom);
			return
	  	end
	end
end

fclose(fecom);

disp(' ## Done.')
