What is the acceptable level of orthogonality errors? I would say none. Or, more
precisely, the question is irrelevant, because orthogonality errors can be kept at
the level roundoff errors of double precision -- it is technologically achievable,
and therefore should be done this way. Then there is no point to think about how
orthogonality errors corrupt the solution.
Regarding mesh: depending on the shape of the domain and coastline curvilinear grids
may be very useful, and are presently underutilized. The main reason, in my view, is
embarrassing lack of necessary software tools to generate such grids and to manipulate
data, pre- and post-processing.
If we ever to build a working sea-grid kind tool, it should be built around Laplace
equation solver with Mehrstellen discretization: it is a 9-point scheme, it is
compact fourth-order accurate. These is no point to have anything less accurate.
Laplace equations for both horizontal coordinates are coupled through boundary
conditions. Currently they are not.
Current/previous/existing approaches try to adapt ready-to-use algorithms, whatever
is available in Matlab or elsewhere. But if nothing fits precisely, then basically
settle for surrogates or whatever means to have it half-way, or simply give up.
An example? bilinear interpolation problem. Suppose point (xc,yc) is somewhere inside
the trapezoidal element defined by 4 points (i,j), (i+1,j), (i+1,j+1) (i,j+1) of the
"parent" (generally speaking orthogonal curvilinear) grid; then there must exist two
number p,q bounded by [0,1] each, such that
Code:
xc=(1-p)*(1-q)*x(i,j) + p*(1-q)*x(i+1,j) + (1-p)*q*x(i,j+1) +p*q*x(i+1,j+1)
yc=(1-p)*(1-q)*y(i,j) + p*(1-q)*y(i+1,j) + (1-p)*q*y(i,j+1) +p*q*y(i+1,j+1)
assuming that xc,yc, and all x(i,j),y(i,j) are known, find p, and q.
How to solve it?
ROMS community answer: it cannot be done.
It is nonlinear problem.
Matlab does not provide it.
Therefore we use griddata - we have no other choice.
Griddata splits the trapezoidal element into two triangles. How does it split? By one
of the diagonals. Which one? Either one, 50% chance (b.t.w, the diagonals are exactly
equal to each other, if it is orthogonal curvilinear grid). Then in each triangle it can
uniquely define a linear function of the coordinates x,y and solve for the coefficients.
This is not exactly what we want, but this is what we settle for.
And before one gets to the above problem of finding p,q one should find indices (i,j)
of the parent grid cell which contains point (xc,yc). This requires search. A not so
trivial problem if the grid is curvilinear.
Again, Matlab does not provide it.
so, griddata.
Quote:
Regarding mesh, I can see that many wants to nest inside global models (Hycom or
NEMO) and do not want/need to download the whole grid domain just to get stripes
along the boundary. In other words, getting stripe along constant lon, and/or lat of global
model saves disk/time using opendap subsettings if your ROMS grid is along that line
So desire to have straight lines along latitude or longitude comes from the fact that it
is easy to cut out data from a global dataset, which which is in lon-lat coordinates.
Getting the whole global data is out of question because of its size.
Frankly I would cut-out a logically-rectangular region which fully contains my ROMS model
domain and use 2D interpolation with rotation of velocity vectors as necessary.
In addition, in terms of read-write performance of netCDF files, getting 1D strip vs. getting
2D patch takes comparable time:
netCDF-3: assuming that longitude corresponds to 1st netCDF dimension, latitude to the 2nd,
reading 1D strip along longitude (1st dimension) takes almost no time: you read from the disk
only what you need; in contrast, getting 1D strip along 2nd dimension is comparable to reading
the entire 2D data bounded by the minimum and maximum latitudes of the desired strip.
netCDF4: assuming the same about 1st and 2nd dimensions: getting 1D strips in lon and lat
direction takes about the same time and requires reading the set of blocks which contain
all the data. Simply put, your reading becomes reading of 2D data regardless along which
direction you need your 1D strip. Depending on the netCDF/HDF5 block size vs. your ROMS
grid geographical extents, it may take exactly the same time as to cut out the entire 2D
region, or less.