seagrid2roms.m

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

nc{'h'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'h'}.long_name = ncchar('Final bathymetry at RHO-points');
nc{'h'}.units = ncchar('meter');
nc{'h'}.field = ncchar('bath, scalar');
 
nc{'f'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'f'}.long_name = ncchar('Coriolis parameter at RHO-points');
nc{'f'}.units = ncchar('second-1');
nc{'f'}.field = ncchar('Coriolis, scalar');
 
nc{'pm'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'pm'}.long_name = ncchar('curvilinear coordinate metric in XI');
nc{'pm'}.units = ncchar('meter-1');
nc{'pm'}.field = ncchar('pm, scalar');
 
nc{'pn'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'pn'}.long_name = ncchar('curvilinear coordinate metric in ETA');
nc{'pn'}.units = ncchar('meter-1');
nc{'pn'}.field = ncchar('pn, scalar');
 
nc{'dndx'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'dndx'}.long_name = ncchar('xi derivative of inverse metric factor pn');
nc{'dndx'}.units = ncchar('meter');
nc{'dndx'}.field = ncchar('dndx, scalar');
 
nc{'dmde'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'dmde'}.long_name = ncchar('eta derivative of inverse metric factor pm');
nc{'dmde'}.units = ncchar('meter');
nc{'dmde'}.field = ncchar('dmde, scalar');
 
nc{'x_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'x_rho'}.long_name = ncchar('x location of RHO-points');
nc{'x_rho'}.units = ncchar('meter');
 
nc{'y_rho'} = ncdouble('eta_rho', 'xi_rho'); %% 16900 elements.
nc{'y_rho'}.long_name = ncchar('y location of RHO-points');
nc{'y_rho'}.units = ncchar('meter');
 
nc{'x_psi'} = ncdouble('eta_psi', 'xi_psi'); %% 16641 elements.
nc{'x_psi'}.long_name = ncchar('x location of PSI-points');
nc{'x_psi'}.units = ncchar('meter');
 
nc{'y_psi'} = ncdouble('eta_psi', 'xi_psi'); %% 16641 elements.
nc{'y_psi'}.long_name = ncchar('y location of PSI-points');
nc{'y_psi'}.units = ncchar('meter');
 
nc{'x_u'} = ncdouble('eta_u', 'xi_u'); %% 16770 elements.
nc{'x_u'}.long_name = ncchar('x location of U-points');
nc{'x_u'}.units = ncchar('meter');
 
nc{'y_u'} = ncdouble('eta_u', 'xi_u'); %% 16770 elements.
nc{'y_u'}.long_name = ncchar('y location of U-points');
nc{'y_u'}.units = ncchar('meter');
 
nc{'x_v'} = ncdouble('eta_v', 'xi_v'); %% 16770 elements.
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');
nc{'angle'}.units = ncchar('radian');   %yes, ROMS expects radians

% Fill the variables with data.

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

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

% Fill the variables.

% Need (x..., y...) in meters.  Currently, they are
%  in Mercator (projected) units.

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

nc{'xl'}(:) = xl;
nc{'el'}(:) = el;

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

nc{'x_rho'}(:) = grid_x(j_rho, i_rho);
nc{'y_rho'}(:) = grid_y(j_rho, i_rho);

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

if ~isempty(bathymetry)
		 		  nc{'h'}(:) = bathymetry;
end

% Masking.

mask = ~~mask;
land = mask;
water = ~land;

rmask = water;

% Calculate other masking arrays.

umask = zeros(size(rmask));
vmask = zeros(size(rmask));
pmask = zeros(size(rmask));

for i = 2:LP
  for j = 1:MP
    umask(j, i-1) = rmask(j, i) * rmask(j, i-1);
  end
end

for i = 1:LP
  for j = 2:MP
    vmask(j-1, i) = rmask(j, i) * rmask(j-1, i);
  end
end

for i = 2:LP
  for j = 2:MP
    pmask(j-1, i-1) = rmask(j, i) * rmask(j, i-1) * rmask(j-1, i) * rmask(j-1, i-1);
  end
end

nc{'mask_rho'}(:) = rmask;
nc{'mask_psi'}(:) = pmask(1:end-1, 1:end-1);
nc{'mask_u'}(:) = umask(1:end, 1:end-1);
nc{'mask_v'}(:) = vmask(1:end-1, 1:end);

% Average angle -- We should do this via (x, y) components.

temp = 0.5*(ang(1:end-1, :) + ang(2:end, :));
ang = zeros(n, m);
ang(2:2:end, 2:2:end) = temp;

nc{'angle'}(:) = ang(j_rho, i_rho);

if (0)
  sx = abs(spaced_x(:, 2:end) - spaced_x(:, 1:end-1));
  sy = abs(spaced_x(2:end, :) - spaced_x(1:end-1, :));
elseif (0)
  sx = abs(spaced_x(1:2:end, 3:2:end) - spaced_x(1:2:end, 1:2:end-2));
  sy = abs(spaced_y(3:2:end, 1:2:end) - spaced_y(1:2:end-2, 1:2:end));
  sx = 0.5 * (sx(1:end-1, :) + sx(2:end, :));
  sy = 0.5 * (sy(:, 1:end-1) + sy(:, 2:end));
end

% Use geometry from seagrid file.
% Note: need half the number of points.

gx = geometry{1};   % Spherical distances in meters.
gy = geometry{2};

% raw_grid_size = [m, n]

% geometry_sizes = [size(gx) size(gy)]

sx = 0.5*(gx(1:end-1, :) + gx(2:end, :));
sy = 0.5*(gy(:, 1:end-1) + gy(:, 2:end));

% raw_s_sizes = [size(sx) size(sy)]

% sx = sx(2:end-1, :);
% sy = sy(:, 2:end-1);

pm = 1 ./ sx;
pn = 1 ./ sy;

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;

% Final size of file:

s = size(nc);
disp([' ## Dimensions: ' int2str(s(1))])
disp([' ## Variables: ' int2str(s(2))])
disp([' ## Global Attributes: ' int2str(s(3))])
disp([' ## Record Dimension: ' name(recdim(nc))])
 
endef(nc)
close(nc)