All,
As we look more closely into re() and compare it with the classic
regrid2() the issue of "polar corrections" came to the surface. As I
understand it, both re() and regrid2() currently apply a polar
correction which consists of replacing the polar values with a single
constant value equal to the average of the interpolated values.
Mike/BJ: please jump in if I am missing something.
Jennifer/Brian: how is this handled by linterp()?
Now, the current polar correction is fine for scalars like temperature
or wind speed that are unique at the poles. This is not true, however,
for the (u,v) wind components since the unit vectors change with
longitude. There is a simple formula for this that can be derived by
writing the wind vector in cartesian coordinates at the poles and
imposing uniqueness of the wind velocity (I have code somewhere for
doing this.) The main point is that you need to have both u & v in
order to apply this correction.
On the same subject, let's digress a bit about polar interpolation...
Case 1: Input fields contain the poles
-----------------------------------------------
Now, when the input field include the poles, there is no latitudinal
interpolation taking place, only interpolation in longitude. Rather
then force u & v to be constant over the poles (this is really wrong,
messes up vorticity/divergence terribly), it is better to apply no
correction. This would work also for scalars: the input field would
have constant values over the pole, so there would be no interpolation
error in longitude.
The upshot: if the input grid contains the poles, apply no polar
correction whatsoever. One can always provide a polar correction
functions:
ufixpole(u,v)
vfixpole(u,v)
if you really would like the poles to be 100% consistent, say to
ensure that the wind speed is unique.
Case 2: input fields do not contain the poles, but output fields do
--------------------------------------------------------------------------------
Now, when the input fields DO NOT contain the poles (say, offset by
half delta-lat), but the *output* fields contain the pole, one should
resist the urge to extrapolate from lower latitudes. We can still do
interpolation because the sphere is doubly periodic after all. A
common technique is to use hemispheric continuation to find the wind
components at the other side of the poles. For example, if last
latitude is
latJ = 90 - dLat/2
then linear interpolation would give us
u(90N,lon) = ( u(latJ,lon) - u(latJ,lon+180) ) / 2
v(90N,lon) = ( v(latJ,lon) - v(latJ,lon+180) ) / 2
(notice the minus sign once you cross the pole - this is needed
because the unit vectors change directions when you are looking from
the other side of the pole.) For a scalar one would have
h(90N,lon) = ( h(latJ,lon) + h(latJ,lon+180) ) / 2.
(notice the plus sign in this case --- the are no unit vectors
involved here.) Similar formulas hold for the South Pole, and for
bessel interpolation with an additional point across the pole being
involved.)
Pole correction: you would need it here. But it may be better to
implement this separately because it is different for vectors and
scalars. The same ufixpole/vfixpole above would be applicable, and a
specific function for scalars.
Proposal
-----------
In retrospect, here is perhaps what would make more sense. We need 2
interpolation functions (or a single function with an optional
parameter selecting the desired behavior -- I lean towards presenting
the user with 2 separate function as it makes it very clear the 2
distinct cases):
re_scalar(expr,...)
re_vector(expr,...)
Now, the scalar version would always apply the pole correction (simply
replace the polar value with the average value as it is currently
done). The re_vector() function would apply no polar correction, and
functions
ufixpole(u,v)
vfixpole(u,v)
would be provided to fix the poles since it needs both u & v to do
this. A similar approach would apply to linterp() as well.
Any comments/suggestions greatly appreciated.
Arlindo
--
Arlindo da Silva
da...@al...
|
