# 1 "fish.c"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "fish.c"
# 312 "fish.c"
=pod

=head1 NAME

fish.gex - GrADS Extension Library for Calculating Streamfunction/Velocity Potential

=head1 SYNOPSIS

=head3 GrADS Functions:

=over 4

=item

B<fish>(I<UEXPR,VEXPR[,MBDCND]>) - Poisson Solver

=item

B<fish_psi>(I<UEXPR,VEXPR[,MBDCND]>) - Computes Stream Function

=item

B<fish_chi>(I<UEXPR,VEXPR[,MBDCND]>0 - Computes Velocity Potential

=item

B<fish_vor>(I<UEXPR,VEXPR[,MBDCND]>) - Computes Relative Vorticity

=item

B<fish_div>(I<UEXPR,VEXPR[,MBDCND]>0 - Computes Divergence

=back

=head1 DESCRIPTION

This library provides GrADS extensions (I<gex>) with functions for
computation of streamfunction and velocity potential from zonal and
meridional wind components:

   laplacian(psi) = vorticity (1)
   laplacian(chi) = divergence (2)

where C<psi> is the streamfunction and C<chi> is the velocity
potential. (See Wikipedia links below for more information on
streamfunction/velocity potential.) The vorticity and divergence
computation relies on functions B<madvu> and B<madvv> provided in the
OpenGrADS extension I<Libbjt> by B.-J. Tsuang. The Poisson equations
(1)-(2) above are solved using the classic C<Fishpak> Fortran
library. Documentation for C<FISHPAK> is given in:

Swarztrauber, P. and R. Sweet, 1975: Efficient Fortran Subprograms for
the Solution of Elliptic Partial Differential Equations. I<Technical Note
TN/IA-109>. National Center for Atmospheric Research, Boulder, Colorado 80307.

All functions provided require a global, uniform lat/lon grid.

=head1 EXAMPLES

For the following examples it is suggested that you open the following
OPeNDAP dataset with a GFS forecast. Fire up C<gradsdods> and then

   ga-> sdfopen http:

=head2 Computing Streamfunction from Vorticity

If you have the relative vorticity field available, say C<vor>, you
can easily compute streamfunction

   ga-> psi = fish(vor)

=head2 Computing Streamfunction from Wind Components

The first step is to evaluate the streamfunction:

  ga-> set lev 200
  ga-> psi = fish_psi(ugrd,vgrd)

It is often convenient to I<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

=head2 Computing the Rotational Wind

Although the intrinsc GrADS function C<cdiff> could be used for
numerically diffferentiating the streamfunction, we use here functions
C<mvadv/muadv> from I<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

=head2 Computing Velocity Potential from Divergence

If you have the divergence field available, say C<div>, you
can easily compute streamfunction

   ga-> chi = fish(div)

=head2 Computing Velocity Potential from Wind Components

Start by computing the velocity potential:

  ga-> set lev 200
  ga-> chi = fish_chi(ugrd,vgrd)

It is often convenient to I<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

=head2 Calculating the Divergent Wind

Although the intrinsc GrADS function C<cdiff> could be used for
numerically diffferentiating the velocity potential, we use here
functions C<mvadv/muadv> from I<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

=head2 Computing Vorticity and Divergence

As an alternative to the intrinsic GrADS functions C<hcurl/hdivg>,
functions C<fish_vor> and C<fish_div> uses the advention functions in
I<Libbjt> to numerically evaluate vorticity and divergence

   ga-> vor = fish_vor(ugrd,vgrd)
   ga-> div = fish_div(ugrd,vgrd)

These functions provide a better handling of the boundaries compared
to their intrinsic counterparts.

=head1 FUNCTIONS PROVIDED

=head2 B<fish>(I<UEXPR,VEXPR[,MBDCND]>) - Poisson Solver

=over 4

This function uses routine C<PWSSSP> in C<Fishpak> to solve the
poisson equation. The default parameter C<MBDCND=9> solves a Helmholts
equation with a very small random term to provide a "unique" solution.

=over 8

=item I<UEXPR,VEXPR> - required

GrADS expressions with zonal and meridional wind components

=item I<MBDCND> - optional

Meridional boundary condition; the default C<MBDCND=9> should work in
most cases. See C<FISHPAK> documentation for additional information.

      MBDCND = 1 ! BC: solution specified at both poles
      MBDCND = 5 ! BC: solution specified at TF (South Pole) and
      MBDCND = 7 ! BC: solution specified at TS (North Pole) and
      MBDCND = 9 ! BC: solution unspecified at both poles

=back

=back

=head2 B<fish_psi>(I<UEXPR,VEXPR[,MBDCND]>) - Computes Stream Function

=over 4

This function computes vorticity as in B<fish_vor> and uses B<fish> to
solve the Poisson equation for the streamfunction C<psi>:

   laplacian(psi) = vorticity

=over 8

=item I<UEXPR,VEXPR> - required

GrADS expressions with zonal and meridional wind components

=item I<MBDCND> - optional

Meridional boundary condition; the default C<MBDCND=9> should work in
most cases. See C<FISHPAK> documentation for additional information.

      MBDCND = 1 ! BC: solution specified at both poles
      MBDCND = 5 ! BC: solution specified at TF (South Pole) and
      MBDCND = 7 ! BC: solution specified at TS (North Pole) and
      MBDCND = 9 ! BC: solution unspecified at both poles

=back

=back

=head2 B<fish_chi>(I<UEXPR,VEXPR[,MBDCND]>) - Computes Velocity Potential

=over 4

This function computes divergence as in B<fish_div> and uses B<fish> to
solve the Poisson equation for the velocity potential C<chi>:

   laplacian(chi) = divergence

=over 8

=item I<UEXPR,VEXPR> - required

GrADS expressions with zonal and meridional wind components

=item I<MBDCND> - optional

Meridional boundary condition; the default C<MBDCND=9> should work in
most cases. See C<FISHPAK> documentation for additional information.

      MBDCND = 1 ! BC: solution specified at both poles
      MBDCND = 5 ! BC: solution specified at TF (South Pole) and
      MBDCND = 7 ! BC: solution specified at TS (North Pole) and
      MBDCND = 9 ! BC: solution unspecified at both poles

=back

=back


=head2 B<fish_vor>(I<UEXPR,VEXPR[,MBDCND]>0 - Computes Relative Vorticity

=over 4

This function computes the vorticity using the expression:

   vorticity = - ( madvu(v,one) - madvv(u,cosphi) / cosphi )

where C<u> and C<v> are the zonal/meridional wind components, C<one>
is a constant field equal to 1, and C<cosphi> is the cosine of
latitude. The functions B<madvu> and B<madvv> are provided in the OpenGrADS
extension library I<Libbjt>.

=over 8

=item I<UEXPR,VEXPR> - required

GrADS expressions with zonal and meridional wind components

=item I<MBDCND> - optional

Meridional boundary condition; the default C<MBDCND=9> should work in
most cases. See C<FISHPAK> documentation for additional information.

      MBDCND = 1 ! BC: solution specified at both poles
      MBDCND = 5 ! BC: solution specified at TF (South Pole) and
      MBDCND = 7 ! BC: solution specified at TS (North Pole) and
      MBDCND = 9 ! BC: solution unspecified at both poles

=back

=back


=head2 B<fish_div>(I<UEXPR,VEXPR[,MBDCND]>0 - Computes Divergence

=over 4

This function computes the divergence using the expression:

   divergence = - ( madvu(u,one) + madvv(v,cosphi) / cosphi )

where C<u> and C<v> are the zonal/meridional wind components, C<one>
is a constant field equal to 1, and C<cosphi> is the cosine of
latitude. The functions B<madvu> and B<madvv> are provided in the OpenGrADS
extension library I<Libbjt>.


=over 8

=item I<UEXPR,VEXPR> - required

GrADS expressions with zonal and meridional wind components

=item I<MBDCND> - optional

Meridional boundary condition; the default C<MBDCND=9> should work in
most cases. See C<FISHPAK> documentation for additional information.

      MBDCND = 1 ! BC: solution specified at both poles
      MBDCND = 5 ! BC: solution specified at TF (South Pole) and
      MBDCND = 7 ! BC: solution specified at TS (North Pole) and
      MBDCND = 9 ! BC: solution unspecified at both poles

=back

=back

=head1 BUGS

The function C<fish> assumes a global and uniform longitude/latitude grid. If
your grid is not uniform, consider using C<re()> for interpolating to
a uniform grid. When doing so, do not make the poles gridpoints. (The
function C<muadv> from I<Libbjt> produces undefined values if the first
and last latitudinal gridpoints are at the poles.)

Undefined values are handled in a less than idea way by the Poisson
solver C<fish>. If undefined values of vortivity/divergence occur at
the first and last latitudinal gridpoint, the following I<polar fix>
is applied depending on whether these gridpoints are the poles or not:

=over 4

=item Poles are gridpoints

If the first and last latitudinal gridpoints are at the poles, the
zonal average of the adjascent latitudinal band is computed and this
zonal averaged value is used at the pole.

=item Poles are not gridpoints

The value at the same longitude in the adjascent latitudinal band is
used.

=back

For interior points, undefined values are set to zero before solving
the Poisson equation.


=head1 SEE ALSO

=over 4

=item *

L<http:

=item *

L<http:

=item *

L<http:

=item *

L<http:

=item *

L<http:
definition on Wikipedia.

=back

=head1 AUTHOR

Arlindo da Silva (dasilva@opengrads.org)

=head1 COPYRIGHT

Copyright (C) 2007-2008 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.

=cut
