c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                          .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c This file contains documentation for subroutine lfp.  The subroutine
c is in the file alf.f
c
c subroutine lfp (init,n,m,l,cp,pb,w)
c
c dimension of           cp(n/2+1), pb(l), w(5*l+41)
c arguments
c
c purpose                routine lfp uses coefficients computed by
c                        routine alfk to calculate the single precision
c                        normalized associated legendre function pbar(n,
c                        m,theta) at colatitudes theta=(i-1)*pi/(l-1),
c                        i=1,...,l. subroutine lfp evaluates pbar
c                        using one of the following trigonometric
c                        expansions
c
c                        1) for n even and m even, pbar(m,n,theta) =
c                           .5*cp(1) plus the sum from k=1 to k=n/2
c                           of cp(k+1)*cos(2*k*th)
c
c                        2) for n even and m odd, pbar(m,n,theta) =
c                           the sum from k=1 to k=n/2 of
c                           cp(k)*sin(2*k*th)
c
c                        3) for n odd and m even, pbar(m,n,theta) =
c                           the sum from k=1 to k=(n+1)/2 of
c                           cp(k)*cos((2*k-1)*th)
c
c                        4) for n odd and m odd,  pbar(m,n,theta) =
c                           the sum from k=1 to k=(n+1)/2 of
c                           cp(k)*sin((2*k-1)*th)
c
c
c usage                  call lfp(init,n,m,l,cp,pb,w)
c
c arguments
c
c on input               init
c                          = 0 initialization only
c                          = 1 compute pbar(n,m,theta)
c
c                          lfp call with init = 0 initializes array w;
c                          no values of pbar(n,m,theta) are computed.
c                          init=0 should be used on the first call, or
c                          if l or w values differ from those in the
c                          previous call.
c
c                        n
c                          nonnegative integer, less than l, specifying
c                          the degree of pbar(n,m,theta)
c
c                        m
c                          is the order of pbar(n,m,theta). m can be
c                          any integer however pbar(n,m,theta) = 0
c                          if abs(m) is greater than n and
c                          pbar(n,m,theta) = (-1)**m*pbar(n,-m,theta)
c                          for negative m.
c
c                        l
c                          number of colatitudes theta=(i-1)*pi/(l-1)
c                          for i=1,...,l where l is greater than 1.
c                          l must be an odd integer.
c
c                        cp
c                          single precision array of length (n/2)+1
c                          containing coefficients computed by routine
c                          alfk
c
c                        w
c                          a single precision work array with at
c                          least 5*l+41 locations
c
c on output              pb
c                          single precision array of length l containing
c                          pbar(n,m,theta), theta=(i-1)*pi/(l-1) for i=1
c                          ,...,l.
c
c                        w
c                          a single precision array containing values
c                          which must not be destroyed if the next call
c                          will have the same value of input parameter n
c
c special conditions     calls to routine lfp must be preceded by an
c                        appropriate call to routine alfk.
c
c precision              single
c
c algorithm              the trigonometric series formula used by
c                        routine lfp to calculate pbar(n,m,theta) for
c                        theta=(i-1)*pi/(l-1), i=1,...,n, depends on
c                        m and n as follows:
c
c                           1) for n even and m even, the formula is
c                              .5*cp(1) plus the sum from k=1 to k=n/2
c                              of cp(k+1)*cos(2*k*theta)
c                           2) for n even and m odd. the formula is
c                              the sum from k=1 to k=n/2 of
c                              cp(k)*sin(2*k*theta)
c                           3) for n odd and m even, the formula is
c                              the sum from k=1 to k=(n+1)/2 of
c                              cp(k)*cos((2*k-1)*theta)
c                           4) for n odd and m odd, the formula is
c                              the sum from k=1 to k=(n+1)/2 of
c                              cp(k)*sin((2*k-1)*theta)
c
c accuracy               comparison between routines lfp and double
c                        precision dlfp on the cray1 indicates greater
c                        accuracy for smaller values of input parameter
c                        n.  agreement to 12 places was obtained for
c                        n=10 and to 11 places for n=100.
c
c timing                 time per call to routine lfp is dependent on
c                        the input parameters l and n.
c