| shfilt.gex - GrADS Extension Library with Spherical Harmonic Utilities |
shfilt.gex - GrADS Extension Library with Spherical Harmonic Utilities
display sh_filt(EXPR,N1[,N2]) - Spherical Harmonic Filter
display sh_power(EXPR) - Spherical Harmonic Power Spectra
display sh_fish(EXPR) - Poisson Solver
display sh_psi(UEXPR,VEXPR) - Computes Stream Function
display sh_chi(UEXPR,VEXPR) - Computes Velocity Potential
display sh_vor(UEXPR,VEXPR) - Computes Relative Vorticity
display sh_div(UEXPR,VEXPR) - Computes Divergence
run power.gs EXPR N - Plot Spherical Harmonic Power Spectra
This library provides GrADS extensions (gex) with functions for spherical harmonic filtering and calculation of spherical harmonic power spectra. The numerical calculations rely on the excellent Spherepak library by John C. Adams and Paul N. Swarztrauber. Quoting from the Speherepak website:
SPHEREPACK 3.1 is a collection of FORTRAN programs that facilitates computer modeling of geophysical processes. The package contains programs for computing certain common differential operators including divergence, vorticity, gradients, and the Laplacian of both scalar and vector functions. Programs are also available for inverting these operators. For example, given divergence and vorticity, the package can be used to compute the velocity components. The Laplacian can also be inverted and therefore the package can be used to solve both the scalar and vector Poisson equations. Its use in model development is demonstrated by a sample program that solves the time-dependent non-linear shallow-water equations. Accurate solutions are obtained via the spectral method that uses both scalar and vector spherical harmonic transforms that are available to the user. The package also contains utility programs for computing the associated Legendre functions, Gauss points and weights, and multiple fast Fourier transforms. Programs are provided for both equally-spaced and Gauss distributed latitudinal points as well as programs that transfer data between these grids.
The current GrADS extensions only begin to explore the capabilties of Spherepak. The function sh_filt takes a scalar global 2D field on the sphere, expands it in terms of spherical harmonics, and reconstructs it includng only (total) wavenumbers in the range [N1,N2] specified by the user. Additionally, function sh_power returns the power spectra in terms of total waenumbers. This is accomplished by returning a 1D array (fixed longitude, varying latitude) with the spectra as a function of total wavenumber. The GrADS script power.gs is useful to plot this power spectra and should be used in place of the function sh_power().
This library also provides functions for computation of streamfunction and velocity potential from zonal and meridional wind components using a spectral method:
laplacian(psi) = vorticity (1) laplacian(chi) = divergence (2)
where psi is the streamfunction and chi is the velocity
potential. (See Wikipedia links below for more information on
streamfunction/velocity potential.) The vorticity and divergence
computation relies on functions madvu and madvv provided in the
OpenGrADS extension Libbjt by B.-J. Tsuang. The Poisson equations
(1)-(2) above are solved using a spectral method in Spherepak. These
functions complement similar functions in the OpenGrADS fish extension
which become numerically unstable for resolutions finer than 1/2 deg.
The example expands the surface pressure field in terms of spherical harmonics and reconstructs it retaining only 10 wavenumber.
open model set gxout shaded d sh_filt(ps,10)
The example expands the surface temperature field in terms of spherical harmonics and plots it as function of the 32 first total numbers.
run power.gs ts 32
If you have the relative vorticity field available, say vor, you
can easily compute streamfunction
ga-> psi = sh_fish(vor)
The first step is to evaluate the streamfunction:
ga-> set lev 200 ga-> psi = sh_psi(ugrd,vgrd)
It is often convenient to center the streamfunction by subtracting its global mean:
ga-> psi = psi - aave(psi,global)
We can finally display it:
ga-> set gxout shaded ga-> display psi/1e7 draw title Streamfunction
Notice that if you would like to compare the results of
sh_psi and fish_psi you must first remove the global
mean of both.
Although the intrinsc GrADS function cdiff could be used for
numerically diffferentiating the streamfunction, we use here functions
mvadv/muadv from Libbjt because of its handling of the
boundaries:
ga-> one = 1 + 0 * lat ga-> upsi = mvadv(one,psi) ga-> vpsi = - muadv(one,psi)
We then plot the rotational streamlines on top of the streamfunction:
ga-> set gxout shaded ga-> display psi/1e7 ga-> set gxout stream ga-> display upsi;vpsi ga-> draw title Streamfunction/Rotational Wind
If you have the divergence field available, say div, you
can easily compute streamfunction
ga-> chi = sh_fish(div)
Start by computing the velocity potential:
ga-> set lev 200 ga-> chi = sh_chi(ugrd,vgrd)
It is often convenient to center the velocity potential by subtracting its global mean:
ga-> psi = psi - aave(psi,global)
We can finally display it:
ga-> set gxout shaded ga-> display chi/1e6 draw title Velocity
Notice that if you would like to compare the results of
sh_chi and fish_chi you must first remove the global
mean of both.
Although the intrinsc GrADS function cdiff could be used for
numerically diffferentiating the velocity potential, we use here
functions mvadv/muadv from Libbjt because of its handling of the
boundaries:
ga-> uchi = - muadv(one,chi) ga-> vchi = - mvadv(one,chi)
We finally display the divergent wind as vectors on top of the velocity potential:
ga-> display chi/1e6 ga-> set gxout vector ga-> set cmin 2 ga-> set cmax 20 ga-> display skip(uchi,6,6);vchi ga-> draw title Velocity Potential/Divergent Wind
As an alternative to the intrinsic GrADS functions hcurl/hdivg,
sh_vor and sh_div uses the advention functions in
Libbjt to numerically evaluate vorticity and divergence
ga-> vor = sh_vor(ugrd,vgrd) ga-> div = sh_div(ugrd,vgrd)
These functions provide a better handling of the boundaries compared to their intrinsic counterparts.
You cannot hcurl/hdivg with sh_fish, you must use
sh_vor/sh_div instead.
This function takes a scalar global 2D field on the sphere, expands it in terms of spherical harmonics, and reconstructs it includng only (total) wavenumbers in the range [N1,N2] specified by the user.
GrADS expressions with scalar expression to be filtered.
When both N1 and N2 are specified, only spherical harmonics with total wavenumber in the range [N1,N2] will be retained. When only N1 is specified, only spherical harmonics with total wavenumber in the range [1,N1] will be retained.
This function returns the power spectra in terms of total waenumbers. This is accomplished by returning a 1D array (fixed longitude, varying latitude) with the spectra as a function of total wavenumber. The GrADS script power.gs is useful to plot this power spectra and should be used in place of the function sh_power().
This function returns the inverse laplacian of EXPR using a spectral method. This is particularly useful when the finite difference implementation in fish() fails for high resolution grids.
This function computes vorticity as in sh_vor and uses sh_fish to
solve the Poisson equation for the streamfunction psi:
laplacian(psi) = vorticity
GrADS expressions with zonal and meridional wind components
This function computes divergence as in sh_div and uses sh_fish to
solve the Poisson equation for the velocity potential chi:
laplacian(chi) = divergence
GrADS expressions with zonal and meridional wind components
This function computes the vorticity using the expression:
vorticity = - ( madvu(v,one) - madvv(u,cosphi) / cosphi )
where u and v are the zonal/meridional wind components, one
is a constant field equal to 1, and cosphi is the cosine of
latitude. The functions madvu and madvv are provided in the OpenGrADS
extension library Libbjt.
GrADS expressions with zonal and meridional wind components
This function computes the divergence using the expression:
divergence = - ( madvu(u,one) + madvv(v,cosphi) / cosphi )
where u and v are the zonal/meridional wind components, one
is a constant field equal to 1, and cosphi is the cosine of
latitude. The functions madvu and madvv are provided in the OpenGrADS
extension library Libbjt.
GrADS expressions with zonal and meridional wind components
Do not use the GrADS intrinsic functions hcurl and hdivg with sh_fish. Use the provided functions sh_vor and sh_div instead.
http://opengrads.org/ - OpenGrADS Home Page
http://cookbooks.opengrads.org/index.php - OpenGrADS Cookbooks
http://www2.cisl.ucar.edu/resources/legacy/spherepack - Spherepack Home Page
http://opengrads.org/wiki/index.php - OpenGrADS User Defined Extensions
http://www.iges.org/grads/ - Official GrADS Home Page
http://en.wikipedia.org/wiki/Velocity_potential - Velocity Potential definition on Wikipedia.
http://en.wikipedia.org/wiki/Stream_function - Stream function definition on Wikipedia.
Arlindo da Silva (dasilva@opengrads.org) with contributions from Ben Auer (Poisson Solver.)
Copyright (C) 2008-2009 Arlindo da Silva; All Rights Reserved.
This is free software; see the source for copying conditions. There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
| shfilt.gex - GrADS Extension Library with Spherical Harmonic Utilities |