seagrid2scrum.m

一些制作正交曲线网格的matlab源程序

function seagrid2scrum(varargin)

% function seagrid2scrum(theSeagridFile, theScrumFile, theGridTitle)

% seagrid2scrum -- Output SCRUM format from SeaGrid file.
%  seagrid2scrum('theSeagridFile', 'theScrumFile', 'theGridTitle')
%   creates a SCRUM file, based on the given SeaGrid file.  The
%   get/put dialog boxes are invoked where filenames are absent
%   or given by a wildcard.
%
% NOTE: this routine does not employ _FillValue attributes
%  in the output NetCDF variables.
 
% 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 18-Jun-1999 17:11:59.
% Updated    31-May-2000 10:33:39.

if (1)
	disp([' ## seagrid2scrum is obsolete; calling "seagrid2roms" instead.'])
	if nargin < 1
		feval('seagrid2roms')
	else
		feval('seagrid2roms', varargin{:})
	end
	return
end

if nargin < 1, theSeagridFile = '*.mat'; end
if nargin < 2, theScrumFile = 'scrum_model_grid.nc'; end
if nargin < 3, theGridTitle = char(zeros(1, 128)+abs(' ')); end

if isempty(theSeagridFile) | 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 | isempty(theSeagridFile) | any(theScrumFile == '*')
	[f, p] = uiputfile(theScrumFile, 'Save to Scrum File:');
	if ~any(f), return, end
	if p(end) ~= filesep, p(end+1) = filesep; end
	theScrumFile = [p f]
end

disp([' ## SeaGrid Source File  : ' theSeagridFile])
disp([' ## Scrum Destination File: ' theScrumFile])

% Load the SeaGrid file and get parameters.

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

% With grid_x of size [m, n], the grid itself has
%  [m-1, n-1] cells.  The latter size corresponds
%  to the size of the mask and bathymetry.  These
%  cell-centers are called the "rho" points.

s = theSeagridData.s;

grid_x = s.grids{1};
grid_y = s.grids{2};
[m, n] = size(grid_x);

geogrid_lon = s.geographic_grids{1};
geogrid_lat = s.geographic_grids{2};
geometry = s.geometry;
mask = s.mask;   % land = 1; water = 0.
if ~isequal(size(mask), size(grid_x)-1)
	if ~isempty(mask)
		disp(' ## Wrong size mask.')
	end
	mask = zeros(m-1, n-1);
end

bathymetry = s.gridded_bathymetry;
projection = s.projection;
ang = s.orientation;
min_depth = s.clipping_depths(1);
max_depth = s.clipping_depths(2);
coriolis = geogrid_lon;

spaced_x = s.spaced_grids{1};
spaced_y = s.spaced_grids{2};

% Double the grid-size before proceeding.
%  The grid-cell-centers are termed the "rho" points.

theInterpFcn = 'interp2';
theInterpMethod = 'spline';

grid_x = feval(theInterpFcn, grid_x, 1, theInterpMethod);
grid_y = feval(theInterpFcn, grid_y, 1, theInterpMethod);
geogrid_lon = feval(theInterpFcn, geogrid_lon, 1, theInterpMethod);
geogrid_lat = feval(theInterpFcn, geogrid_lat, 1, theInterpMethod);
spaced_x = feval(theInterpFcn, spaced_x, 1, theInterpMethod);
spaced_y = feval(theInterpFcn, spaced_y, 1, theInterpMethod);

% The present size of the grid nodes.

[m, n] = size(grid_x);

% Flip arrays top for bottom.
%  Is this necessary?

FLIPPING = 0;
FLIPPING = 1;

if FLIPPING
	grid_x = flipud(grid_x);
	grid_y = flipud(grid_y);
	geogrid_lon = flipud(geogrid_lon);
	geogrid_lat = flipud(geogrid_lat);
	geometry{1} = flipud(geometry{1});
	mask = flipud(mask);
	bathymetry = flipud(bathymetry);
	ang = flipud(ang);
	coriolis = flipud(coriolis);
	spaced_x = flipud(spaced_x);
	spaced_y = flipud(spaced_y);
end

% Create the Scrum NetCDF file.

%% ncdump('priapus:MATLAB:toolbox:scrum:tasman_grd.nc')   %% Generated 17-Jun-1999 17:37:30
 
nc = netcdf(theScrumFile, 'clobber');
if isempty(nc), return, end
 
%% Global attributes:

disp(' ## Defining Global Attributes...')
 
nc.type = ncchar('Gridpak file');
nc.gridid = theGridTitle;
nc.history = ncchar(['Created by "' mfilename '" on ' datestr(now)]);

nc.CPP_options = ncchar('DCOMPLEX, DBLEPREC, NCARG_32, PLOTS,');
name(nc.CPP_options, 'CPP-options')

% The SeaGrid is now a full array, whose height
%  and width are odd-valued.  We extract staggered
%  sub-grids for the SCRUM scheme, ignoring the
%  outermost rows and columns.  Thus, the so-called
%  "rho" points correspond to the even-numbered points
%  in an (i, j) Matlab array.  The "psi" points begin
%  at i = 3 and j = 3.  The whole set is indexed as
%  follows:

% rho (2:2:end-1, 2:2:end-1), i.e. (2:2:m, 2:2:n), etc.
% psi (3:2:end-2, 3:2:end-2)
% u   (2:2:end-1, 3:2:end-2)
% v   (3:2:end-2, 2:2:end-1)

if ~rem(m, 2), m = m-1; end   % m, n must be odd.
if ~rem(n, 2), n = n-1; end

i_rho = 2:2:m-1; j_rho = 2:2:n-1;
i_psi = 3:2:m-2; j_psi = 3:2:n-2;
i_u = 2:2:m-1; j_u = 3:2:n-2;
i_v = 3:2:m-2; j_v = 2:2:n-1;

% The xi direction (left-right):
%	LP = (m-1)/2;   % The rho dimension.
LP = (n-1)/2;   % The rho dimension.
L = LP-1;       % The psi dimension.

% The eta direction (up-down):
%	MP = (n-1)/2;   % The rho dimension.
MP = (m-1)/2;   % The rho dimension.
M = MP-1;       % The psi dimension.

[L M LP MP]

disp(' ## Defining Dimensions...')
 
nc('xi_psi') = L;
nc('xi_rho') = LP;
nc('xi_u') = L;
nc('xi_v') = LP;

nc('eta_psi') = M;
nc('eta_rho') = MP;
nc('eta_u') = MP;
nc('eta_v') = M;

nc('two') = 2;
nc('bath') = 0; %% (record dimension)
 
%% Variables and attributes:

disp(' ## Defining Variables and Attributes...')
 
nc{'xl'} = ncdouble; %% 1 element.
nc{'xl'}.long_name = ncchar('domain length in the XI-direction');
nc{'xl'}.units = ncchar('meter');
 
nc{'el'} = ncdouble; %% 1 element.
nc{'el'}.long_name = ncchar('domain length in the ETA-direction');
nc{'el'}.units = ncchar('meter');
 
nc{'JPRJ'} = ncchar('two'); %% 2 elements.
nc{'JPRJ'}.long_name = ncchar('Map projection type');

nc{'JPRJ'}.option_ME_ = ncchar('Mercator');
nc{'JPRJ'}.option_ST_ = ncchar('Stereographic');
nc{'JPRJ'}.option_LC_ = ncchar('Lambert conformal conic');
name(nc{'JPRJ'}.option_ME_, 'option(ME)')
name(nc{'JPRJ'}.option_ST_, 'option(ST)')
name(nc{'JPRJ'}.option_LC_, 'option(LC)')
 
nc{'PLAT'} = ncfloat('two'); %% 2 elements.
nc{'PLAT'}.long_name = ncchar('Reference latitude(s) for map projection');
nc{'PLAT'}.units = ncchar('degree_north');
 
nc{'PLONG'} = ncfloat; %% 1 element.
nc{'PLONG'}.long_name = ncchar('Reference longitude for map projection');
nc{'PLONG'}.units = ncchar('degree_east');
 
nc{'ROTA'} = ncfloat; %% 1 element.
nc{'ROTA'}.long_name = ncchar('Rotation angle for map projection');
nc{'ROTA'}.units = ncchar('degree');
 
nc{'JLTS'} = ncchar('two'); %% 2 elements.
nc{'JLTS'}.long_name = ncchar('How limits of map are chosen');

nc{'JLTS'}.option_CO_ = ncchar('P1, .. P4 define two opposite corners ');
nc{'JLTS'}.option_MA_ = ncchar('Maximum (whole world)');
nc{'JLTS'}.option_AN_ = ncchar('Angles - P1..P4 define angles to edge of domain');
nc{'JLTS'}.option_LI_ = ncchar('Limits - P1..P4 define limits in u,v space');
name(nc{'JLTS'}.option_CO_, 'option(CO)')
name(nc{'JLTS'}.option_MA_, 'option(MA)')
name(nc{'JLTS'}.option_AN_, 'option(AN)')
name(nc{'JLTS'}.option_LI_, 'option(LI)')
 
nc{'P1'} = ncfloat; %% 1 element.
nc{'P1'}.long_name = ncchar('Map limit parameter number 1');
 
nc{'P2'} = ncfloat; %% 1 element.
nc{'P2'}.long_name = ncchar('Map limit parameter number 2');
 
nc{'P3'} = ncfloat; %% 1 element.
nc{'P3'}.long_name = ncchar('Map limit parameter number 3');
 
nc{'P4'} = ncfloat; %% 1 element.
nc{'P4'}.long_name = ncchar('Map limit parameter number 4');
 
nc{'XOFF'} = ncfloat; %% 1 element.
nc{'XOFF'}.long_name = ncchar('Offset in x direction');
nc{'XOFF'}.units = ncchar('meter');
 
nc{'YOFF'} = ncfloat; %% 1 element.
nc{'YOFF'}.long_name = ncchar('Offset in y direction');
nc{'YOFF'}.units = ncchar('meter');
 
nc{'depthmin'} = ncshort; %% 1 element.
nc{'depthmin'}.long_name = ncchar('Shallow bathymetry clipping depth');
nc{'depthmin'}.units = ncchar('meter');
 
nc{'depthmax'} = ncshort; %% 1 element.
nc{'depthmax'}.long_name = ncchar('Deep bathymetry clipping depth');
nc{'depthmax'}.units = ncchar('meter');
 
nc{'spherical'} = ncchar; %% 1 element.
nc{'spherical'}.long_name = ncchar('Grid type logical switch');
nc{'spherical'}.option_T_ = ncchar('spherical');
nc{'spherical'}.option_F_ = ncchar('Cartesian');
name(nc{'spherical'}.option_T_, 'option(T)')