support for reduced gaussian grid/glesp pixelization

[Former user]

Hi. I've been using shtns to develop a weather prediction model

gfs-dycore.googlecode.com

It's a version of the model that is used by the U.S. National Weather Service. The operational NWS version uses a grid called in the weather community a "reduced gaussian grid", which I believe is the astronomy community is called "Gauss- Legendre Sky Pixelization" (GLESP). The idea behind reduced gaussian grids is described here

http://journals.ametsoc.org/doi/abs/10.1175/1520-0493%281991%29119%3C1057:UORGGI%3E2.0.CO%3B2

and the GLESP pixelization is described here

http://www.glesp.nbi.dk/

The basic idea is to use the Gaussian-Legendre quadrature points in the latitude direction, but with a variable number of points in the longitude direction (decreasing from 2*NLATS at the equator toward zero at the poles so that the grid is approximately equal area).

I was wondering if you had any plans to include this kind of grid as an option in SHTNS. It would result in a computational savings of 20-30% in the global weather model case. I would be willing to work on this, if you could give me a rough outline of what would need to be done.

Thanks for you work on SHTNS, it's a great library and a terrific resource for global weather modelling.

Regards, Jeff

[Nathanaël Schaeffer]

(Reply via [e-mail]):

Hello,

Well, the reduced grids you propose will effectively speed up your calculations only if you spend a lot of time computing things in the physical domain. Is that the case for your application ? Because if not, polar-optimization is implemented in SHTns, which neglects those high m spherical harmonics at high latitudes.

I would guess that a reduced grid will not make the transform faster, because currently the fft's are treated together efficiently by FFTW, and if you end up with different number of grid points in longitude, you'll have to treat them independently, which will probably lead to longer execution time.

So if you convince me that it is interesting (because you spend so much time computing stuff in the physical domain), I will consider to implement it . Otherwise, a rather quick implementation would only require you to change the Fourier transform part (by splitting it up), which may not be so easy because I pack to latitudinal rings into one complex number series to get faster Fourier transform. You could then in addition adjust the array shtns->tm[im] which contains the index to the first latidude to consider for each m.

Please, tell me if you are receiving this email, since I'm not yet used to this issue-tracker.

Regards, Nat.

[Former user]

Nate: Weather models spend a lot of time computing things like the effects of clouds, turbulence and radiation in the physical domain - that's way cutting down the number of grid points can speed things up. If the reduced gaussian grid will not speed up the transforms (or will in fact slow them down) I could just implement an extra interpolation step to go back and forth between the full and reduced grids everytime a transform is done. -Jeff

[Nathanaël Schaeffer]

Ok, so are there some constraints on the way the spatial data should be stored, or am I free to do anything I want ? (I would of course provide a map giving cos(theta) and phi for every data point) In fact, I'd like to store theta and pi-theta together. Does that sound ok ?

[Former user]

Storing as a 1d array would be fine.