[Opengrads-devel] Fishy finite differencing

[Arlindo da S. <da...@al...>]

Mike,

  Or should I have said "finite differencing for fish()"?. (This is intended
primary for Mike; I'm cc'ing the list in case someone wants to comment on
the  algorithm.)

   A little background. Fish() is a poisson solver that computes
streamfunction and velocity potential given vorticity and divervence:

   psi = fish(vort)
   chi = fish(divg)

Unlike the classic UDF psichi,  I  have  separated  the computation of
vort/div from the possion solver proper. Like the classic psichi,  fish()
uses fishpak.f, and old dusty-deck f77 eliptic solver that stood the test of
time. (Fishpak is not really f77 but rather f66, with tons of goto's and
arithmetic ifs, but it is very fast and the numerics are solid.) Fishpak is
freely available from netlib.

Now, in principle vorticity and divergence could be computed with grads
intrinsics such as:

  vort = hcurl(u,v)

The problem is, hcurl() does not handles the boundaries very well, placing
UNDEFS in the first and last points (even if this is a global periodic
grid). Altough fish replaces undefs with 0. (like the classic psichi), this
is not the right thing to do. So, we need to provide a replacement in the
fish() package for properly handling the boundaries, something that we could
contribute later to COLA for improving the current intrinsics.

Here is my proposal for handling finite differencing at the boundaries:

a) Zonal direction. If the grid is global, exploit the periodicidy, and use
centered differences for all points. If the grid (read: dimension
environment) is regional, use uncentered differences at the boundaries.

b) Meridional direction. If the grid (dimension environment) is regional,
use uncentered differences at the boundaries. Now, if the grid is global,
the handling of the extrema depends on whether the pole is a grid point or
not:

i) Grid is off the pole by delta-lat/2. For example, on a 1 degree grid the
first and last gridpoints are -89.5 and 89.5. In this case use uncentered
differences at these points.

ii) The poles are grid-points. This is the grid used by the finite-volume
dynamical core of S.J.-Lin used at GMAO, NCAR and GFDL.  The trick here is
to consider these polar points as little polar caps and use vector calculus
to compute vorticity/divergence. Using Stokes' theorem one can relate the
mean vorticity in this polar cap to the "circulation" at the boundaries of
the polar cap, basically the sum of "u"  at the boundary of the polar cap
for some 2 pi factor thrown in for good measure; for 1 deg grids, this
boundary is at -89.5 and 89.5, so "u" needs to be interpolated as the
average of the first and second meridion grid points (similarly for last and
1 from last gridpoints).  A similar algorithm is used for divergence, but
now one relies on the Gauss (divergence) theorem to relate the mean
divergence at the polar cap to the mass flux into the polar cap, basically
the sum of "v" at the boundary with some 2pi normalization factor.

iii) Gaussian grids. The best solution would be to do a spectral transform
and compute vort/div that way.

Poles are always messy, and usually a big source of trouble. So it pays to
be a bit careful.

Mike: I added the stubs already to the fish() code base, perhaps more than
we need. I guess all that we need at this point is vort() and div(). Could
you add your code for the zonal direction? You can implement the polar
polar, and if you wish I can do it or write down the algorithm more
precisely (or send you the references.)

I need to dig some spectral transform code for the gaussian grid case, or
else do a simple uncentered calculation  at the pole. If one is going to
bring spectral transforms in, then it is better to do the calculation of
psi/chi spectrally and not use fish() which is intended for lat/lon grids.
This may be another extension altogether.

I'm cc'ing Brian in case he has any suggestions for us, as he may want to
incorporate some of this into hcur/divg later.

    As always, any suggestions greatly appreciated.

      Cheers!

         Arlindo
















-- 
Arlindo da Silva
da...@al...