mat2roms.m

任意边界结构正交曲线网格生成程序

nc{'x_v'}.long_name = ncchar('x location of V-points');
nc{'x_v'}.units = ncchar('meter');
 
nc{'y_v'} = ncdouble('eta_v', 'xi_v'); %% 16770 elements.
nc{'y_v'}.long_name = ncchar('y location of V-points');
nc{'y_v'}.units = ncchar('meter');
 
nc{'lat_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'lat_rho'}.long_name = ncchar('latitude of RHO-points');
nc{'lat_rho'}.units = ncchar('degree_north');
 
nc{'lon_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'lon_rho'}.long_name = ncchar('longitude of RHO-points');
nc{'lon_rho'}.units = ncchar('degree_east');
 
nc{'lat_psi'} = ncdouble('eta_psi', 'xi_psi'); %% 16641 elements.
nc{'lat_psi'}.long_name = ncchar('latitude of PSI-points');
nc{'lat_psi'}.units = ncchar('degree_north');
 
nc{'lon_psi'} = ncdouble('eta_psi', 'xi_psi'); %% 16641 elements.
nc{'lon_psi'}.long_name = ncchar('longitude of PSI-points');
nc{'lon_psi'}.units = ncchar('degree_east');
 
nc{'lat_u'} = ncdouble('eta_u', 'xi_u'); %% 16770 elements.
nc{'lat_u'}.long_name = ncchar('latitude of U-points');
nc{'lat_u'}.units = ncchar('degree_north');
 
nc{'lon_u'} = ncdouble('eta_u', 'xi_u'); %% 16770 elements.
nc{'lon_u'}.long_name = ncchar('longitude of U-points');
nc{'lon_u'}.units = ncchar('degree_east');
 
nc{'lat_v'} = ncdouble('eta_v', 'xi_v'); %% 16770 elements.
nc{'lat_v'}.long_name = ncchar('latitude of V-points');
nc{'lat_v'}.units = ncchar('degree_north');
 
nc{'lon_v'} = ncdouble('eta_v', 'xi_v'); %% 16770 elements.
nc{'lon_v'}.long_name = ncchar('longitude of V-points');
nc{'lon_v'}.units = ncchar('degree_east');
 
nc{'mask_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'mask_rho'}.long_name = ncchar('mask on RHO-points');
nc{'mask_rho'}.option_0_ = ncchar('land');
nc{'mask_rho'}.option_1_ = ncchar('water');
name(nc{'mask_rho'}.option_0_, 'option(0)')
name(nc{'mask_rho'}.option_1_, 'option(1)')
 
nc{'mask_u'} = ncdouble('eta_u', 'xi_u'); %% 16770 elements.
nc{'mask_u'}.long_name = ncchar('mask on U-points');
nc{'mask_u'}.option_0_ = ncchar('land');
nc{'mask_u'}.option_1_ = ncchar('water');
name(nc{'mask_u'}.option_0_, 'option(0)')
name(nc{'mask_u'}.option_1_, 'option(1)')
%		 nc{'mask_u'}.FillValue_ = ncdouble(1);
 
nc{'mask_v'} = ncdouble('eta_v', 'xi_v'); %% 16770 elements.
nc{'mask_v'}.long_name = ncchar('mask on V-points');
nc{'mask_v'}.option_0_ = ncchar('land');
nc{'mask_v'}.option_1_ = ncchar('water');
name(nc{'mask_v'}.option_0_, 'option(0)')
name(nc{'mask_v'}.option_1_, 'option(1)')
%		 nc{'mask_v'}.FillValue_ = ncdouble(1);
 
nc{'mask_psi'} = ncdouble('eta_psi', 'xi_psi'); %% 16641 elements.
nc{'mask_psi'}.long_name = ncchar('mask on PSI-points');
nc{'mask_psi'}.option_0_ = ncchar('land');
nc{'mask_psi'}.option_1_ = ncchar('water');
name(nc{'mask_psi'}.option_0_, 'option(0)')
name(nc{'mask_psi'}.option_1_, 'option(1)')
%		 nc{'mask_psi'}.FillValue_ = ncdouble(1);

% Now, what about depths: "h" and "hraw".  <== DEPTHS.

% The following seems mistaken:
 
nc{'angle'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'angle'}.long_name = ncchar('angle between xi axis and east');
nc{'angle'}.units = ncchar('degree');

% Variables.

% Fill the variables with data.

disp(' ## Filling Variables...')

if ~isfield(s, 'projection') | isempty(s.projection)
	s.projection = 'Mercator';
end
projection = s.projection;

switch lower(projection)
case 'mercator'
	theProjection = 'ME';
case 'stereographic'
	theProjection = 'ST';
case 'lambert conformal conic'
	theProjection = 'LC';
otherwise
	theProjection = '??';
end

nc{'JPRJ'}(1:2) = theProjection;
nc{'spherical'}(:) = 'T';   % T or F -- uppercase okay?

% Cartesian coordinates.

x = s.rho.x;
y = s.rho.y;
grid_x = interp2(x, 1);
grid_y = interp2(y, 1);

% Geographic coordinates.

meters = 1/(40e6);   % Earth circumference (reciprocal).

m_proj(s.projection)
[lon, lat] = m_xy2ll(x, y);
[geogrid_lon, geogrid_lat] = m_xy2ll(grid_x*meters, grid_y*meters);   % Degrees.

xl = max(grid_x(:)) - min(grid_x(:));
el = max(grid_y(:)) - min(grid_y(:));
nc{'xl'}(:) = xl;
nc{'el'}(:) = el;

f = 2.*7.29e-5.*sin(lat.*pi./180);
nc{'f'}(:) = f;

% Handy indices.

[m, n] = size(grid_x);

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

% Locations.

nc{'x_rho'}(:) = s.rho.x;
nc{'y_rho'}(:) = s.rho.y;

nc{'x_psi'}(:) = grid_x(j_psi, i_psi);
nc{'y_psi'}(:) = grid_y(j_psi, i_psi);

nc{'x_u'}(:) = grid_x(j_u, i_u);
nc{'y_u'}(:) = grid_y(j_u, i_u);

nc{'x_v'}(:) = grid_x(j_v, i_v);
nc{'y_v'}(:) = grid_y(j_v, i_v);

nc{'lon_rho'}(:) = geogrid_lon(j_rho, i_rho);
nc{'lat_rho'}(:) = geogrid_lat(j_rho, i_rho);

nc{'lon_psi'}(:) = geogrid_lon(j_psi, i_psi);
nc{'lat_psi'}(:) = geogrid_lat(j_psi, i_psi);

nc{'lon_u'}(:) = geogrid_lon(j_u, i_u);
nc{'lat_u'}(:) = geogrid_lat(j_u, i_u);

nc{'lon_v'}(:) = geogrid_lon(j_v, i_v);
nc{'lat_v'}(:) = geogrid_lat(j_v, i_v);

% Metric factors: use ones if absent.

if ~isfield(s.rho, 'dx') | isempty(s.rho.dx)
	s.rho.dx = ones(size(nc{'pm'}));
end
if ~isfield(s.rho, 'dy') | isempty(s.rho.dy)
	s.rho.dy = ones(size(nc{'pn'}));
end

dx = s.rho.dx;
dy = s.rho.dy;

dx(dx == 0) = NaN;   % Shouldn't be any zeros here.
dy(dy == 0) = NaN;

pm = 1 ./ dx;   % Note reciprocals.
pn = 1 ./ dy;

% Does ROMS tolerate NaNs and/or Infs?
%  Do we have to substitute a NetCDF "_FillValue"?

nc{'pm'}(:) = pm;
nc{'pn'}(:) = pn;

dmde = zeros(size(pm));
dndx = zeros(size(pn));

dmde(2:end-1, :) = 0.5*(1./pm(3:end, :) - 1./pm(1:end-2, :));
dndx(:, 2:end-1) = 0.5*(1./pn(:, 3:end) - 1./pn(:, 1:end-2));

nc{'dmde'}(:) = dmde;
nc{'dndx'}(:) = dndx;

% Grid-cell angles, in degrees: compute if necessary.

if ~isfield(s.rho, 'angle') | isempty(s.rho.angle)
	dx = diff(s.rho.x, 1, 2);
	dx(1, :) = [];
	dy = -diff(s.rho.y, 1, 1);
	dy(:, end) = [];
	ang = angle(dx+dy*sqrt(-1)) * 180 / pi;
	ang = [ang(1, :); ang];   % Pad at top.
	ang = [ang ang(:, end)];   % Pad on right.
	s.rho.angle = ang;
end

nc{'angle'}(:) = s.rho.angle;   % Degrees.

% Depths at RHO points.

bathymetry = s.rho.depth;
mask = isnan(bathymetry) | (bathymetry == -99999);   % ECOM land = -99999.
bathymetry(mask) = 0;
bathymetry(mask) = NaN;   % Which is better?  ROMS: NaNs okay?
water = 1 - mask;   % ROMS requires water = 1, land= 0.
if ~isempty(bathymetry)
	nc{'h'}(:) = bathymetry;   % ROMS depths are positive.
end

% The u, v, and psi masks rely directly on
%  the mask of the surrounding rho points.

nc{'mask_rho'}(:) = water;
nc{'mask_u'}(:) = water(:, 1:end-1) & water(:, 2:end);
nc{'mask_v'}(:) = water(1:end-1, :) & water(2:end, :);
nc{'mask_psi'}(:) = water(1:end-1, 1:end-1) &...
									water(1:end-1, 2:end) & ...
									water(2:end, 1:end-1) & ...
									water(2:end, 2:end);

% Close the ROMS File.

if ~isempty(close(nc))
	disp(' ## Unable too close the ROMS output file.')
end

% WetCDF off.

if isMacintosh
	eval('wetcdf off')
end

% ---------- romsnan ---------- %

function r = romsnan

% romsnan -- NaN value prefered by ROMS.
%  romsnan (no argument) returns the value prefered
%   by ROMS in place of NaN.  *** I believe it is
%   the NetCDF _FillValue for doubles.  (CRD)

f = ncfillvalues;
r = f.ncdouble;