Abstract

Library FFTPACK 5.1 contains 1D, 2D, and multiple fast Fourier subroutines, written in Fortran 77, for transforming real and complex data, real even and odd wave data, and real even and odd quarter-wave data. All of the FFTPACK 5.1 routines listed above are grouped in triplets e.g. {CFFT1I, CFFT1F, CFFT1B}. The suffix I denotes initialize, F denotes forward (as in forward transform) and B denotes backward. In an application program, before calling B or F routines for the first time, or before calling them with a different length, users must initialize an array by calling the I routine of the appropriate pair or triplet. Note that I routines need not be called each time before a B or F routine is called.

All of the transform routines in FFTPACK 5.1 are normalized.

Error messages are written to unit 6 by routine XERFFT. The standard version of XERFFT issues an error message and halts execution, so that no FFTPACK 5.1 routine will return to the calling program with error return IER different than zero. Users may consider modifying the STOP statement in order to call system-specific exception-handling facilities.

Caveat

FFTPACK 5.1 is not fully compliant with the Fortran 77 standard. There are several instances where arrays of type REAL or COMPLEX are passed to a subroutine and used as a different type. We have not passed the code through a rigorous standards-checker, so we do not have a list of the infractions. Prospective users whose projects require strict adherence to the Fortran standard should not use this library.

References

(1) Vectorizing the Fast Fourier Transforms, by Paul Swarztrauber, Parallel Computations, G. Rodrigue, ed., Academic Press, New York 1982.

(2) Fast Fourier Transforms Algorithms for Vector Computers, by Paul Swarztrauber, Parallel Computing, (1984) pp.45-63.

Overview

Complex Transform Routines
__________________________

CFFT1I    1D complex initialization
CFFT1B    1D complex backward
CFFT1F    1D complex forward

CFFT2I    2D complex initialization
CFFT2B    2D complex backward
CFFT2F    2D complex forward

CFFTMI    multiple complex initialization
CFFTMB    multiple complex backward
CFFTMF    multiple complex forward


Real Transform Routines
_______________________

RFFT1I    1D real initialization
RFFT1B    1D real backward
RFFT1F    1D real forward 

RFFT2I    2D real initialization
RFFT2B    2D real backward
RFFT2F    2D real forward 

RFFTMI    multiple real initialization
RFFTMB    multiple real backward
RFFTMF    multiple real forward 


Real Cosine Transform Routines
______________________________

COST1I    1D real cosine initialization
COST1B    1D real cosine backward
COST1F    1D real cosine forward 

COSTMI    multiple real cosine initialization
COSTMB    multiple real cosine backward
COSTMF    multiple real cosine forward 


Real Sine Transform Routines
____________________________

SINT1I    1D real sine initialization
SINT1B    1D real sine backward
SINT1F    1D real sine forward 

SINTMI    multiple real sine initialization
SINTMB    multiple real sine backward
SINTMF    multiple real sine forward 


Real Quarter-Cosine Transform Routines
______________________________________

COSQ1I    1D real quarter-cosine initialization
COSQ1B    1D real quarter-cosine backward
COSQ1F    1D real quarter-cosine forward 

COSQMI    multiple real quarter-cosine initialization
COSQMB    multiple real quarter-cosine backward
COSQMF    multiple real quarter-cosine forward 


Real Quarter-Sine Transform Routines
____________________________________

SINQ1I    1D real quarter-sine initialization
SINQ1B    1D real quarter-sine backward
SINQ1F    1D real quarter-sine forward 

SINQMI    multiple real quarter-sine initialization
SINQMB    multiple real quarter-sine backward
SINQMF    multiple real quarter-sine forward 

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFT1I (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)

DESCRIPTION

 FFTPACK 5.1 subroutine CFFT1I initializes array WSAVE for use in 
 its companion routines CFFT1B and CFFT1F.  Routine CFFT1I must 
 be called before the first call to  CFFT1B or CFFT1F, and after 
 whenever the value of integer N changes.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG(REAL(N))/LOG(2.)) + 4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines CFFT1B or CFFT1F.


 IER     = 0 successful exit
         = 2 input parameter LENSAV not big enough

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFT1B (N, INC, C, LENC, WSAVE, LENSAV,
1                  WORK, LENWRK, IER)

 INTEGER    N, INC, LENC, LENSAV, LENWRK, IER
 COMPLEX    C(LENC)
 REAL       WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine CFFT1B computes the one-dimensional Fourier 
 transform of a single periodic sequence within a complex array.  
 This transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequence from spectral to physical 
 space.
 
 This transform is normalized since a call to CFFT1B followed
 by a call to CFFT1F (or vice-versa) reproduces the original
 array subject to algorithm constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array C, of two 
         consecutive elements within the sequence to be transformed.
 
 C       Complex array of length LENC containing the sequence to be 
         transformed.
 
 LENC    Integer dimension of C array.  LENC must be at least 
         INC*(N-1) + 1.


 WSAVE   Real work array with dimension LENSAV.  WSAVE's contents 
         must be initialized with a call to subroutine CFFT1I before 
         the first call to routine CFFT1F or CFFT1B for a given 
         transform length N.  WSAVE's contents may be re-used for 
         subsequent calls to CFFT1F and CFFT1B with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG(REAL(N))/LOG(2.)) + 4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at
         least 2*N.


 Output Arguments
 
 C       For index J*INC+1 where J=0,...,N-1,
 
            C(J*INC+1) = 
 
            N-1
            SUM C(K*INC+1)*EXP(I*J*K*2*PI/N)
            K=0

         where I=SQRT(-1).
 
         At other indices, the output value of C does not differ
         from input.
 
 IER     =  0 successful exit
         =  1 input parameter LENC   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFT1F (N, INC, C, LENC, WSAVE, LENSAV,
1                  WORK, LENWRK, IER)

 INTEGER    N, INC, LENC, LENSAV, LENWRK, IER
 COMPLEX    C(LENC)
 REAL       WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine CFFT1F computes the one-dimensional Fourier 
 transform of a single periodic sequence within a complex array.  
 This transform is referred to as the forward transform or Fourier 
 analysis, transforming the sequence from physical to spectral 
 space.
 
 This transform is normalized since a call to CFFT1F followed
 by a call to CFFT1B (or vice-versa) reproduces the original
 array subject to algorithm constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array C, of two 
         consecutive elements within the sequence to be transformed.
 
 C       Complex array of length LENC containing the sequence to be 
         transformed.
 
 LENC    Integer dimension of C array.  LENC must be at least 
         INC*(N-1) + 1.

 WSAVE   Real work array with dimension LENSAV.  WSAVE's contents 
         must be initialized with a call to subroutine CFFT1I before 
         the first call to routine CFFT1F or CFFT1B for a given 
         transform length N.  WSAVE's contents may be re-used for 
         subsequent calls to CFFT1F and CFFT1B with the same N.

 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG(REAL(N))/LOG(2.)) + 4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at
         least 2*N.


 Output Arguments
 
 C       For index J*INC+1 where J=0,...,N-1 (that is, for the Jth 
         element of the sequence),
 
            C(J*INC+1) = 
 
            N-1
            SUM C(K*INC+1)*EXP(-I*J*K*2*PI/N)
            K=0

         where I=SQRT(-1).

         At other indices, the output value of C does not differ
         from input.
 
 IER     =  0 successful exit
         =  1 input parameter LENC   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFT2I (L, M, WSAVE, LENSAV, IER)

 INTEGER    L, M, LENSAV, IER
 REAL       WSAVE(LENSAV)

DESCRIPTION

 FFTPACK 5.1 routine CFFT2I initializes real array WSAVE for use
 in its companion routines CFFT2F and CFFT2B for computing two-
 dimensional fast Fourier transforms of complex data.  Prime 
 factorizations of L and M, together with tabulations of the 
 trigonometric functions, are computed and stored in array WSAVE.
 
 Input Arguments
 
 L       Integer number of elements to be transformed in the first 
         dimension.  The transform is most efficient when L is a 
         product of small primes.


 M       Integer number of elements to be transformed in the second 
         dimension.  The transform is most efficient when M is a 
         product of small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of L and M, and also containing certain 
         trigonometric values which will be used in routines 
         CFFT2B or CFFT2F.


 WSAVE   Real work array with dimension LENSAV.  The WSAVE array 
         must be initialized with a call to subroutine CFFT2I before 
         the first call to CFFT2B or CFFT2F, and thereafter whenever
         the values of L, M or the contents of array WSAVE change.  
         Using different WSAVE arrays for different transform lengths
         or types in the same program may reduce computation costs 
         because the array contents can be re-used.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFT2B (LDIM, L, M, C, WSAVE, LENSAV,
1                     WORK, LENWRK, IER)

 INTEGER    L, M, LDIM, LENSAV, LENWRK, IER
 COMPLEX    C(LDIM,M)
 REAL       WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine CFFT2B computes the two-dimensional discrete
 Fourier transform of a complex periodic array.  This transform is
 known as the backward transform or Fourier synthesis, transforming
 from spectral to physical space.

 Routine CFFT2B is normalized, in that a call to CFFT2B followed
 by a call to CFFT2F (or vice-versa) reproduces the original array
 subject to algorithm constraints, roundoff error, etc.
 
 Input Arguments
 
 LDIM    Integer first dimension of two-dimensional complex array C.
 

 L       Integer number of elements to be transformed in the first 
         dimension of the two-dimensional complex array C.  The value 
         of L must be less than or equal to that of LDIM.  The 
         transform is most efficient when L is a product of small 
         primes.


 M       Integer number of elements to be transformed in the second 
         dimension of the two-dimensional complex array C.  The 
         transform is most efficient when M is a product of small 
         primes.
 
 C       Complex array of two dimensions containing the (L,M) subarray
         to be transformed.  C's first dimension is LDIM, its second 
         dimension must be at least M.


 WSAVE   Real work array with dimension LENSAV.  WSAVE's contents 
         must be initialized with a call to subroutine CFFT2I before 
         the first call to routine CFFT2F or CFFT2B with transform 
         lengths L and M.  WSAVE's contents may be re-used for 
         subsequent calls to CFFT2F and CFFT2B with the same 
         transform lengths L and M.

 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8.

 WORK    Real work array.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at least
         2*L*M.
 
 Output Arguments


  C      Complex output array.  For purposes of exposition,
         assume the index ranges of array C are defined by
         C(0:L-1,0:M-1).
 
         For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given
          in the traditional aliased form by
 
                   L-1  M-1
          C(I,J) = SUM  SUM  C(L1,M1)*
                   L1=0 M1=0
 
                     EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M))
 
          And in unaliased form, the C(I,J)'s are given by
 
                   LF    MF
          C(I,J) = SUM   SUM   C(L1,M1,K1)*
                   L1=LS M1=MS
 
                     EXP(SQRT(-1)*2*PI*(I*L1/L + J*M1/M))
 
          where
 
             LS= -L/2    and LF=L/2-1   if L is even;
             LS=-(L-1)/2 and LF=(L-1)/2 if L is odd;
             MS= -M/2    and MF=M/2-1   if M is even;
             MS=-(M-1)/2 and MF=(M-1)/2 if M is odd;
 
          and
 
             C(L1,M1) = C(L1+L,M1) if L1 is zero or negative;
             C(L1,M1) = C(L1,M1+M) if M1 is zero or negative;
 
          The two forms give different results when used to
          interpolate between elements of the sequence.
 
  IER     Integer error return
          =  0 successful exit
          =  2 input parameter LENSAV not big enough
          =  3 input parameter LENWRK not big enough
          =  5 input parameter L > LDIM
          = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFT2F (LDIM, L, M, C, WSAVE, LENSAV,
1                     WORK, LENWRK, IER)

 INTEGER    L, M, LDIM, LENSAV, LENWRK, IER
 COMPLEX    C(LDIM,M)
 REAL       WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine CFFT2F computes the two-dimensional discrete
 Fourier transform of a complex periodic array.  This transform is
 known as the forward transform or Fourier analysis, transforming
 from physical to spectral space.

 Routine CFFT2F is normalized, in that a call to CFFT2F followed
 by a call to CFFT2B (or vice-versa) reproduces the original array
 subject to algorithm constraints, roundoff error, etc.


 Input Arguments


 LDIM    Integer first dimension of two-dimensional complex array C.


 L       Integer number of elements to be transformed in the first 
         dimension of the two-dimensional complex array C.  The value 
         of L must be less than or equal to that of LDIM.  The 
         transform is most efficient when L is a product of small 
         primes.


 M       Integer number of elements to be transformed in the second 
         dimension of the two-dimensional complex array C.  The 
         transform is most efficient when M is a product of small 
         primes.
 
 C       Complex array of two dimensions containing the (L,M) subarray
         to be transformed.  C's first dimension is LDIM, its second 
         dimension must be at least M.

 WSAVE   Real work array with dimension LENSAV.  WSAVE's contents 
         must be initialized with a call to subroutine CFFT2I before 
         the first call to routine CFFT2F or CFFT2B with transform 
         lengths L and M.  WSAVE's contents may be re-used for 
         subsequent calls to CFFT2F and CFFT2B having those same 
         transform lengths.
 

 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*(L+M) + INT(LOG(REAL(L))/LOG(2.)) + INT(LOG(REAL(M))/LOG(2.)) + 8.


 WORK    Real work array.


 LENWRK  Integer dimension of WORK array.  LENWRK must be at least
         2*L*M.


 Output Arguments


  C      Complex output array.  For purposes of exposition,
         assume the index ranges of array C are defined by
         C(0:L-1,0:M-1).

         For I=0,...,L-1 and J=0,...,M-1, the C(I,J)'s are given
         in the traditional aliased form by

                            L-1  M-1
         C(I,J) =   1/(L*M)*SUM  SUM  C(L1,M1)*
                            L1=0 M1=0

                    EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M))

         And in unaliased form, the C(I,J)'s are given by

                           LF    MF
         C(I,J) =  1/(L*M)*SUM   SUM   C(L1,M1)*
                           L1=LS M1=MS

                    EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M))

         where

            LS= -L/2    and LF=L/2-1   if L is even;
            LS=-(L-1)/2 and LF=(L-1)/2 if L is odd;
            MS= -M/2    and MF=M/2-1   if M is even;
            MS=-(M-1)/2 and MF=(M-1)/2 if M is odd;

         and

            C(L1,M1) = C(L1+L,M1) if L1 is zero or negative;
            C(L1,M1) = C(L1,M1+M) if M1 is zero or negative;

         The two forms give different results when used to
         interpolate between elements of the sequence.

 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  5 input parameter L > LDIM
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFTMI (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)

DESCRIPTION

 FFTPACK 5.1 subroutine CFFTMI initializes array WSAVE for use in 
 its companion routines CFFTMB and CFFTMF.  Routine CFFTMI must 
 be called before the first call to  CFFTMB or CFFTMF, and after 
 whenever the value of integer N changes.
 
 Input Arguments
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG(REAL(N))/LOG(2.)) + 4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines CFFTMB or CFFTMF.


 IER     = 0 successful exit
         = 2 input parameter LENSAV not big enough

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFTMB (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV,
1                  WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER
 COMPLEX    C(LENC)
 REAL       WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine CFFTMB computes the one-dimensional Fourier 
 transform of multiple periodic sequences within a complex array.  
 This transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequences from spectral to physical 
 space.
 
 This transform is normalized since a call to CFFTMF followed
 by a call to CFFTMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array C.
 
 JUMP    Integer increment between the locations, in array C,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array C,
         of two consecutive elements within the same sequence
         to be transformed.
 
 C       Complex array containing LOT sequences, each having 
         length N, to be transformed.  C can have any number
         of dimensions, but the total number of locations must 
         be at least LENC.
 
 LENC    Integer dimension of C array.  LENC must be at 
         least (LOT-1)*JUMP + INC*(N-1) + 1.

 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine CFFTMI before the 
         first call to routine CFFTMF or CFFTMB for a given transform
         length N.  

 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG(REAL(N))/LOG(2.)) + 4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at
         least 2*LOT*N.


 Output Arguments
 
 C       For index L*JUMP+J*INC+1 where J=0,...,N-1 and 
         L=0,...,LOT-1, (that is, for the Jth element of the Lth 
         sequence),
 
            C(L*JUMP+J*INC+1) = 
 
            N-1
            SUM C(L*JUMP+K*INC+1)*EXP(I*J*K*2*PI/N)
            K=0

         where I=SQRT(-1).
 
         At other indices, the output value of C does not differ
         from input.
 
 IER     = 0 successful exit
         = 1 input parameter LENC   not big enough
         = 2 input parameter LENSAV not big enough
         = 3 input parameter LENWRK not big enough
         = 4 input parameters INC,JUMP,N,LOT are not consistent.

             The parameters integers INC, JUMP, N and LOT are 
             consistent if equality 
             I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N 
             and J1,J2 < LOT implies I1=I2 and J1=J2.

             For multiple FFTs to execute correctly, input variables 
             INC, JUMP, N and LOT must be consistent ... otherwise at 
             least one array element mistakenly is transformed more 
             than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE CFFTMF (LOT, JUMP, N, INC, C, LENC, WSAVE, LENSAV,
1                  WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENC, LENSAV, LENWRK, IER
 COMPLEX    C(LENC)
 REAL       WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine CFFTMF computes the one-dimensional Fourier 
 transform of multiple periodic sequences within a complex array.  
 This transform is referred to as the forward transform or Fourier 
 analysis, transforming the sequences from physical to spectral 
 space.
 
 This transform is normalized since a call to CFFTMF followed
 by a call to CFFTMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array C.
 
 JUMP    Integer increment between the locations, in array C,
         of the first elements of two consecutive sequences
         to be transformed.

 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array C,
         of two consecutive elements within the same sequence
         to be transformed.
 
 C       Complex array containing LOT sequences, each having 
         length N, to be transformed.  C can have any number
         of dimensions, but the total number of locations must 
         be at least LENC.
 
 LENC    Integer dimension of C array.  LENC must be at 
         least (LOT-1)*JUMP + INC*(N-1) + 1.

 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine CFFTMI before the 
         first call to routine CFFTMF or CFFTMB for a given transform
         length N. 

 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG(REAL(N))/LOG(2.)) + 4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at
         least 2*LOT*N.


 Output Arguments
 
 C       For index L*JUMP + J*INC +1 where J=0,...,N-1 and 
         L=0,...,LOT-1, (that is, for the Jth element of the Lth 
         sequence),
 
            C(L*JUMP+J*INC+1) = 
 
            N-1
            SUM C(L*JUMP+K*INC+1)*EXP(-I*J*K*2*PI/N)
            K=0

         where I=SQRT(-1).
 
         At other indices, the output value of C does not differ
         from input.
 
 IER     = 0 successful exit
         = 1 input parameter LENC   not big enough
         = 2 input parameter LENSAV not big enough
         = 3 input parameter LENWRK not big enough
         = 4 input parameters INC,JUMP,N,LOT are not consistent.

             The parameters integers INC, JUMP, N and LOT are 
             consistent if equality 
             I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N 
             and J1,J2 < LOT implies I1=I2 and J1=J2.

             For multiple FFTs to execute correctly, input variables 
             INC, JUMP, N and LOT must be consistent ... otherwise at 
             least one array element mistakenly is transformed more 
             than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFT1I (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)

DESCRIPTION

 FFTPACK 5.1 subroutine RFFT1I initializes array WSAVE for use
 in its companion routines RFFT1B and RFFT1F.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines RFFT1B or RFFT1F.


 IER     =  0 successful exit
         =  2 input parameter LENSAV not big enough

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFT1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine RFFT1B computes the one-dimensional Fourier 
 transform of a periodic sequence within a real array.  This  
 is referred to as the backward transform or Fourier synthesis, 
 transforming the sequence from spectral to physical space.
 

 This transform is normalized since a call to RFFT1B followed
 by a call to RFFT1F (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least 
         INC*(N-1) + 1.

 WSAVE   Real work array o length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine RFFT1I before the 
         first call to routine RFFT1F or RFFT1B for a given transform
         length N. 


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(N-1)*INC).
 
         The output values of R are written over the input values.
         If N is even, set NH=N/2-1; then for J=0,...,N-1
 
          R(J*INC) =     R(0) +
 
             [(-1)**J*R((N-1)*INC)]
 
                     NH
               +     SUM  R((2*N1-1)*INC)*COS(J*N1*2*PI/N)
                     N1=1
 
                     NH
               +     SUM  R(2*N1*INC)*SIN(J*N1*2*PI/N)
                     N1=1
 
          If N is odd, set NH=(N-1)/2 and define R as above,
          except remove the expression in square brackets [].
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine RFFT1F computes the one-dimensional Fourier 
 transform of a periodic sequence within a real array.  This  
 is referred to as the forward transform or Fourier analysis, 
 transforming the sequence from physical to spectral space.

 
 This transform is normalized since a call to RFFT1F followed
 by a call to RFFT1B (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least 
         INC*(N-1) + 1.

 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine RFFT1I before the 
         first call to routine RFFT1F or RFFT1B for a given transform
         length N. 


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(N-1)*INC).
 
         Then
 
                        N-1
         R(0) =         SUM  R(N1*INC)/N
                        N1=0
 
         If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2;
         then for J=1,...,NH
 
           R((2*J-1)*INC) = 
 
                     N-1
                  2.*SUM  (R(N1*INC)*COS(J*N1*2*PI/N)/N
                    N1=0
 
          and
 
            R(2*J*INC) = 
 
                     N-1
                  2.*SUM  (R(N1*INC)*SIN(J*N1*2*PI/N)/N
                    N1=0
 
          Also if N is even then
 
            R((N-1)*INC) = 
 
                     N-1
                     SUM  (-1)**N1*R(N1*INC)/N
                     N1=0
 
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFT2I (L, M, WSAVE, LENSAV, IER)
 INTEGER    L, M, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 routine RFFT2I initializes real array WSAVE for use 
 in its companion routines RFFT2F and RFFT2B for computing the two-
 dimensional fast Fourier transform of real data.  Prime 
 factorizations of L and M, together with tabulations of the 
 trigonometric functions, are computed and stored in array WSAVE.
 RFFT2I must be called prior to the first call to RFFT2F or RFFT2B.
 Separate WSAVE arrays are required for different values of L or M.
 
 Input Arguments
 
 L       Integer number of elements to be transformed in the first 
         dimension.  The transform is most efficient when L is a 
         product of small primes.
 
 M       Integer number of elements to be transformed in the second 
         dimension.  The transform is most efficient when M is a 
         product of small primes.
 
 LENSAV  Integer number of elements in the WSAVE array.  LENSAV must
         be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) +
         2*INT(LOG(REAL(M))/LOG(2.)) +12.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of L and M, and also containing certain 
         trigonometric values which will be used in routines 
         RFFT2B or RFFT2F.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFT2B (LDIM, L, M, R, WSAVE, LENSAV, WORK, LENWRK, IER)
 INTEGER    LDIM, L, M, LENSAV, LENWRK, IER
 REAL       R(LDIM,M), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine RFFT2B computes the two-dimensional discrete
 Fourier transform of the complex Fourier coefficients a real 
 periodic array.  This transform is known as the backward transform 
 or Fourier synthesis, transforming from spectral to physical space.
 
 Routine RFFT2B is normalized: a call to RFFT2B followed by a
 call to RFFT2F (or vice-versa) reproduces the original array
 subject to algorithmic sonstraints, roundoff error, etc.
 
 Input Arguments
 
 LDIM    Integer first dimension of the two-dimensional real
         array R, which must be at least L.
 
 L       Integer number of elements to be transformed in the first 
         dimension of the two-dimensional real array R.  The value 
         of L must be less than or equal to that of LDIM.  The 
         transform is most efficient when L is a product of small 
         primes.
 
 M       Integer number of elements to be transformed in the second 
         dimension of the two-dimensional real array R.  The 
         transform is most efficient and accurate when M is a product 
         of small primes.


 R       A real L by M array containing the spectral coefficients
         of a real L by M array that are stored as described in the
         documentation of subroutine rfft2f. THe first dimension is
         LDIM which must be at least L. The second dimension must be 
         at least M.
 

 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine RFFT2I before the 
         first call to routine RFFT2F or RFFT2B. WSAVE's contents may 
         be re-used for subsequent calls to RFFT2F and RFFT2B with 
         the same L and M.


 LENSAV  Integer number of elements in the WSAVE array.  LENSAV must
         be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 
         2*INT(LOG(REAL(M))/LOG(2.)) +12.


 WORK    Real array of dimension LENWRK, where LENWRK is defined 
         below. WORK provides workspace, and its contents need not 
         be saved between calls to routines RFFT2B and RFFT2F.
 
 LENWRK  Integer number of elements in the WORK array.  LENWRK must
         be at least (L+1)*M.

 Output Arguments
 
 R       A real L by M array. If the full transform c is reconstructed 
         using subroutine r2c(ldim,lcdim,l,m,r,c) then for i=0,...,l-1
         and j=0,...,m-1 
            

                      L-1  M-1
             R(I,J) = SUM  SUM  C(L1,M1)
                      L1=0 M1=0

                             *EXP(SQRT(-1)*2*PI*(I*L1/L+J*M1/M))


             or using the conjugate symmetry c(i,j)=c(l-1,m-j) this can
             be written in terms of c(i,j), i=0,...,(L+1)/2 
             and j=0,...,m as:



                            M-1
             R(I,J) = REAL[ SUM C(0,M1)*EXP(SQRT(-1)*2*PI*J*M1/M ]
                            M1=0


                            (L+1)/2-1   M-1
                    + 2*REAL[ SUM       SUM  C(L1,M1)*
                              L1=1      M1=0

                            *EXP(SQRT(-1)*2*PI*(I*L1/L+J*M1/M)) ]

             If L is even then add

                       M-1
             + REAL[ SUM (-1)**I*C(L/2,M1)*EXP(SQRT(-1)*2*PI*J*M1/M) ]
                       M1=0

             c(i,j) = a(i,j)+i*b(I,J) are contained in the real output 
             array r(i,j) except for c(0,m-j) and c(l,m-j) which can
             be obtained from c(0,m-j) = c(0,j) and 
             c(l,j) = c(l.m-j) = c(l,j)



 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  6 input parameter LDIM is less than 2*INT((L+1)/2)
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFT2F (LDIM, L, M, R, WSAVE, LENSAV, WORK, LENWRK, IER)
 INTEGER    LDIM, L, M,  LENSAV, LENWRK, IER
 REAL       R(LDIM,M), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine RFFT2F computes the two-dimensional discrete
 Fourier transform of a real periodic array.  This transform is
 known as the forward transform or Fourier analysis, transforming
 from physical to spectral space.
 
 Routine RFFT2F is normalized: a call to RFFT2F followed by a
 call to RFFT2B (or vice-versa) reproduces the original array
 subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LDIM    Integer first dimension of the two-dimensional real
         array R, which must be at least L.
 
 L       Integer number of elements to be transformed in the first 
         dimension of the two-dimensional real array R.  The value 
         of L must be less than or equal to that of LDIM.  The 
         transform is most efficient when L is a product of small 
         primes.
 
 M       Integer number of elements to be transformed in the second 
         dimension of the two-dimensional real array R.  The 
         transform is most efficient when M is a product of small 
         primes.
 
 R       Real array of two dimensions containing the L-by-M subarray
         to be transformed.  R's first dimension is LDIM and its 
         second dimension must be at least as large as M.
 


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine RFFT2I before the 
         first call to routine RFFT2F or RFFT2B. WSAVE's contents 
         may be re-used for subsequent calls to RFFT2F and RFFT2B 
         as long as L and M remain unchanged.


 LENSAV  Integer number of elements in the WSAVE array.  LENSAV must
         be at least L + 3*M + INT(LOG(REAL(L))/LOG(2.)) + 
         2*INT(LOG(REAL(M))/LOG(2.)) +12.


 WORK    Real array of dimension LENWRK which is defined below.
         WORK provides workspace, and its contents need not be saved 
         between calls to routines RFFT2F and RFFT2B.
 
 LENWRK  Integer number of elements in the WORK array.  LENWRK must
         be at least (L+1)*M.

 Output Arguments
 
 R       Real output array of two dimensions. The full complex transform 
         of r(i,j) is given by:
 
                                 L-1  M-1
              C(I,J) =   1/(L*M)*SUM  SUM  R(L1,M1)*
                                 L1=0 M1=0
   
                    EXP(-SQRT(-1)*2*PI*(I*L1/L + J*M1/M))

         The complex transform of a real array has conjugate symmetry.
         That is, c(i,j) = congugate c(l-i,m-j) so only half the transform
         is computed and packed back into the original array R.

         Examples:   Let a(i,j) = re[c(i,j)] and b(i,j) = Im[c(i,j)] then
         following the forward transform

         For l=m=4
                               a(0,0) a(0,1) b(0,1) a(0,2)
                r(i,j)    =    a(1,0) a(1,1) a(1,2) a(1,3)
                               b(1,0) b(1,1) b(1,2) b(1,3)
                               a(2,0) a(2,1) b(2,1) a(2,2)

         For l=m=5
                               a(0,0) a(0,1) b(0,1) a(0,2),b(0,2)
                               a(1,0) a(1,1) a(1,2) a(1,3),a(1,4)
                r(i,j)     =   b(1,0) b(1,1) b(1,2) b(1,3),b(1,4)
                               a(2,0) a(2,1) a(2,2) a(2,3),a(2,4)
                               b(2,0) b(2,1) b(2,2) b(2,3),b(2,4)

         The remaining c(i,j) for i=int((L+1)/2)+1,..,L and m=0,...m-1
         can be obtained via the conjugate symmetry, which also implies
         that c(0,j) = conjugate c(0,m-j) and for even l, 
         c(l/2,0) = conjugate c(l/2,m-j).

         The full complex transform c(i,j), i=1,...,L and j=1,...,M can also 
         be extracted using
 
                            subroutine r2c(ldim,lcdim,l,m,r,c)

         where lcdim is the first dimension of the complex array c, which 
         must be greater than or equal to l.
         

 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  6 input parameter LDIM is less than 2*INT((L+1)/2)
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFTMI (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)

DESCRIPTION

 FFTPACK 5.1 subroutine RFFTMI initializes array WSAVE for use
 in its companion routines RFFTMB and RFFTMF.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines RFFTMB or RFFTMF.


 IER     =  0 successful exit
         =  2 input parameter LENSAV not big enough

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV,
1                  WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV)     ,WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine RFFTMB computes the one-dimensional Fourier 
 transform of multiple periodic sequences within a real array.  
 This transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequences from spectral to physical 
 space.
 
 This transform is normalized since a call to RFFTMB followed
 by a call to RFFTMF (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at 
         least (LOT-1)*JUMP + INC*(N-1) + 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine RFFTMI before the 
         first call to routine RFFTMF or RFFTMB for a given transform
         length N.  


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at
         least LOT*N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         The output values of R are written over the input values.
         If N is even, set NH=N/2-1; then for I=0,...,LOT-1 and
         J=0,...,N-1
 
          R(I*JUMP+J*INC) =     R(I*JUMP) +
 
             [(-1)**J*R(I*JUMP+(N-1)*INC)]
 
                     NH
               +     SUM  R(I*JUMP+(2*N1-1)*INC)*COS(J*N1*2*PI/N)
                     N1=1
 
                     NH
               +     SUM  R(I*JUMP+2*N1*INC)*SIN(J*N1*2*PI/N)
                     N1=1
 
          If N is odd, set NH=(N-1)/2 and define R as above,
          except remove the expression in square brackets [].
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.

              The parameters integers INC, JUMP, N and LOT are
              consistent if equality
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables
              INC, JUMP, N and LOT must be consistent ... otherwise at
              least one array element mistakenly is transformed more
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE RFFTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV,
1                  WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV)     ,WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine RFFTMF computes the one-dimensional Fourier 
 transform of multiple periodic sequences within a real array.  
 This transform is referred to as the forward transform or Fourier 
 analysis, transforming the sequences from physical to spectral
 space.
 
 This transform is normalized since a call to RFFTMF followed
 by a call to RFFTMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at 
         least (LOT-1)*JUMP + INC*(N-1) + 1.


 WSAVE   Real work array o length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine RFFTMI before the 
         first call to routine RFFTMF or RFFTMB for a given transform
         length N.  
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at
         least LOT*N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         Then for I=0,...,LOT-1
 
                        N-1
         R(I*JUMP) =    SUM  R(I*JUMP+N1*INC)/N
                        N1=0
 
         If N is even, set NH=N/2-1; if N is odd set NH=(N-1)/2;
         then for J=1,...,NH
 
           R(I*JUMP+(2*J-1)*INC) = 
 
                     N-1
                  2.*SUM  (R(I*JUMP+N1*INC)*COS(J*N1*2*PI/N)/N
                    N1=0
 
          and
 
            R(I*JUMP+2*J*INC) = 
 
                     N-1
                  2.*SUM  (R(I*JUMP+N1*INC)*SIN(J*N1*2*PI/N)/N
                    N1=0
 
          Also if N is even then
 
            R(I*JUMP+(N-1)*INC) = 
 
                     N-1
                     SUM  (-1)**N1*R(I*JUMP+N1*INC)/N
                     N1=0
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.

              The parameters integers INC, JUMP, N and LOT are
              consistent if equality
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables
              INC, JUMP, N and LOT must be consistent ... otherwise at
              least one array element mistakenly is transformed more
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COST1I (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine COST1I initializes array WSAVE for use
 in its companion routines COST1F and COST1B.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N-1 is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines COST1B or COST1F.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COST1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COST1B computes the one-dimensional Fourier 
 transform of an even sequence within a real array.  This 
 transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequence from spectral to physical 
 space.
 
 This transform is normalized since a call to COST1B followed
 by a call to COST1F (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N-1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine COST1I before the 
         first call to routine COST1F or COST1B for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to COST1F and COST1B with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least N-1.


 Output Arguments
 
 R       Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(N-1)*INC).
 
         The output values of R are written over the input values.
         For J=0,...,N-1
 
          R(J*INC) =
 
              N-1
              SUM  R(N1*INC)*COS(J*N1*PI/(N-1))
              N1=0
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

SUBROUTINE COST1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)
INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK) 

DESCRIPTION

FFTPACK 5.1 routine COST1F computes the one-dimensional Fourier 
transform of an even sequence within a real array.  This transform 
is referred to as the forward transform or Fourier analysis, 
transforming the sequence from  physical to spectral space. 

This transform is normalized since a call to COST1F followed by a 
call to COST1B (or vice-versa) reproduces the original array 
subject to algorithmic constraints, roundoff error, etc.

Input Arguments 

N       Integer length of the sequence to be transformed.  The transform 
        is most efficient when N-1 is a product of small primes.

INC     Integer increment between the locations, in array R of two consecutive 
        elements within the sequence.

R       Real array of length LENR containing the sequence to be transformed.

LENR    Integer dimension of R array.  LENR must be at least INC*(N-1)+ 1.

WSAVE   Real work array of length LENSAV.  WSAVE's contents must be initialized 
        with a call to subroutine COST1I before the first call to routine COST1F 
        or COST1B for a given transform length N.  WSAVE's contents may be re-used 
        for subsequent calls to COST1F and COST1B with the same N.

LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
        2*N + INT(LOG (REAL(N))/LOG(2.)) +4.

WORK    Real work array of dimension at least LENWRK.

LENWRK  Integer dimension of WORK array.  LENWRK must be at least N-1.

Output Arguments

R       Real output array R.  For purposes of this exposition, assume R's range 
        of indices is given by R(0:(N-1)*INC) and that the input values of R 
        are denoted by X. The output values of R are written over the input 
        values X as follows:

        CASE N=1
        --------
        
        R(0) = X(0)

        CASE N>1
        --------
        
        For J=0

        R(0) = 
             
           0.5*X(0)/(N-1)

           N-2
         + SUM  R(N1*INC)/(N-1)
           N1=1

         + 0.5*X((N-1)*INC)/(N-1)

        For J=1,...,N-2

         R(J*INC) = 

            R(0)/(N-1)

            N-2
          + SUM  2.0*(X(N1*INC)*COS(J*N1*PI/(N-1)))/(N-1)
            N1=1

          + ((-1)**J)*X((N-1)*INC)/(N-1)

         R((N-1)*INC) = 

             0.5*X(0)/(N-1)

             N-2
           + SUM  R(N1*INC)*((-1)**N1)/(N-1)
             N1=1

           + 0.5*((-1)**(N-1))*X((N-1)*INC)/(N-1)

IER     Integer error return
        =  0 successful exit
        =  1 input parameter LENR   not big enough
        =  2 input parameter LENSAV not big enough
        =  3 input parameter LENWRK not big enough
        = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSTMI (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine COSTMI initializes array WSAVE for use
 in its companion routines COSTMF and COSTMB.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines COSTMB or COSTMF.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COSTMB computes the one-dimensional Fourier 
 transform of multiple even sequences within a real array.  This 
 transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequences from spectral to physical 
 space.
 
 This transform is normalized since a call to COSTMB followed
 by a call to COSTMF (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N-1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine COSTMI before the 
         first call to routine COSTMF or COSTMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to COSTMF and COSTMB with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*(N+1).


 Output Arguments
 
 R       Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=0,...,N-1
 
          R(I*JUMP+J*INC) =
 
              N-1
              SUM  R(I*JUMP+N1*INC)*COS(J*N1*PI/(N-1))
              N1=0
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are 
              consistent if equality 
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N 
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables 
              INC, JUMP, N and LOT must be consistent, otherwise at 
              least one array element mistakenly is transformed more 
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COSTMF computes the one-dimensional Fourier 
 transform of multiple even sequences within a real array.  This 
 transform is referred to as the forward transform or Fourier 
 analysis, transforming the sequences from physical to spectral
 space.
 
 This transform is normalized since a call to COSTMF followed
 by a call to COSTMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N-1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine COSTMI before the 
         first call to routine COSTMF or COSTMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to COSTMF and COSTMB with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*(N+1).


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         The output values of R are written over the input values.

         For I=0,...,LOT-1
          R(I*JUMP) = 

              0.5*X(I*JUMP)/(N-1)
 
              N-2
            + SUM  R(I*JUMP+*N1*INC)/(N-1)
              N1=1

            + 0.5*X(I*JUMP+(N-1)*INC)/(N-1)
 
         For I=0,...,LOT-1 and J=1,...,N-2
          R(I*JUMP+J*INC) = 

              R(I*JUMP)/(N-1)
 
              N-2
            + SUM  2.0*(X(I*JUMP+*N1*INC)*COS(J*N1*PI/(N-1)))/(N-1)
              N1=1

            + ((-1)**J)*X(I*JUMP+(N-1)*INC)/(N-1)
 
         For I=0,...,LOT-1
          R(I*JUMP+(N-1)*INC) = 

              0.5*X(I*JUMP)/(N-1)
 
              N-2
            + SUM  R(I*JUMP+*N1*INC)*((-1)**N1)/(N-1)
              N1=1

            + 0.5*((-1)**(N-1))*X(I*JUMP+(N-1)*INC)/(N-1)
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are 
              consistent if equality 
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N 
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables 
              INC, JUMP, N and LOT must be consistent, otherwise at 
              least one array element mistakenly is transformed more 
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINT1I (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine SINT1I initializes array WSAVE for use
 in its companion routines SINT1F and SINT1B.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N+1 is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines SINT1B or SINT1F.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINT1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINT1B computes the one-dimensional Fourier 
 transform of an odd sequence within a real array.  This transform 
 is referred to as the backward transform or Fourier synthesis, 
 transforming the sequence from spectral to physical space.
 
 This transform is normalized since a call to SINT1B followed
 by a call to SINT1F (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N+1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINT1I before the 
         first call to routine SINT1F or SINT1B for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINT1F and SINT1B with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  Must be at least 2*N+2.
 
 Output Arguments
 
  R      Real output array.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(INC:N*INC).
 
         The output values of R are written over the input values.
         For J=1,...,N
 
          R(J*INC) =
 
               N 
              SUM  R(N1*INC)*SIN(J*N1*PI/(N+1))
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINT1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINT1F computes the one-dimensional Fourier 
 transform of an odd sequence within a real array.  This transform 
 is referred to as the forward transform or Fourier analysis, 
 transforming the sequence from physical to spectral space.
 
 This transform is normalized since a call to SINT1F followed
 by a call to SINT1B (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N+1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINT1I before the 
         first call to routine SINT1F or SINT1B for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINT1F and SINT1B with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  Must be at least 2*N+2.
 
 Output Arguments
 
 R       Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by R(INC:(N-1)*INC).
 
         The output values of R are written over the input values.
         For J=1,...,N
 
          R(J*INC) =
 
               N 
              SUM  2.*R(N1*INC)*SIN(J*N1*PI/(N+1))/(N+1)
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINTMI (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine SINTMI initializes array WSAVE for use
 in its companion routines SINTMF and SINTMB.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines SINTMB or SINTMF.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINTMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINTMB computes the one-dimensional Fourier 
 transform of multiple odd sequences within a real array.  This 
 transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequences from spectral to physical 
 space.
 
 This transform is normalized since a call to SINTMB followed
 by a call to SINTMF (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N+1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINTMI before the 
         first call to routine SINTMF or SINTMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINTMF and SINTMB with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*(2*N+4).
 
 Output Arguments
 
  R      Real output array.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(INC:(LOT-1)*JUMP+N*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=1,...,N
 
          R(I*JUMP+J*INC) =
 
              N 
              SUM  R(I*JUMP+*N1*INC)*SIN(J*N1*PI/(N+1))
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are
              consistent if equality
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables
              INC, JUMP, N and LOT must be consistent ... otherwise at
              least one array element mistakenly is transformed more
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINTMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINTMF computes the one-dimensional Fourier 
 transform of multiple odd sequences within a real array.  
 This transform is referred to as the forward transform or Fourier 
 analysis, transforming the sequences from physical to spectral
 space.
 
 This transform is normalized since a call to SINTMF followed
 by a call to SINTMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N+1 is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINTMI before the 
         first call to routine SINTMF or SINTMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINTMF and SINTMB with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         N/2 + N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*(2*N+4).
 
 Output Arguments
 
 R       Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=1,...,N
 
          R(I*JUMP+J*INC) =
 
              N 
              SUM  2.*R(I*JUMP+*N1*INC)*SIN(J*N1*PI/(N+1))/(N+1)
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.

              The parameters integers INC, JUMP, N and LOT are
              consistent if equality
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables
              INC, JUMP, N and LOT must be consistent ... otherwise at
              least one array element mistakenly is transformed more
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSQ1I (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine COSQ1I initializes array WSAVE for use
 in its companion routines COSQ1F and COSQ1B.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines COSQ1B or COSQ1F.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSQ1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COSQ1B computes the one-dimensional Fourier 
 transform of a sequence which is a cosine series with odd wave 
 numbers.  This transform is referred to as the backward transform 
 or Fourier synthesis, transforming the sequence from spectral to 
 physical space.
 
 This transform is normalized since a call to COSQ1B followed
 by a call to COSQ1F (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer number of elements to be transformed in the 
         sequence. The transform is most efficient when N is a 
         product of small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine COSQ1I before the 
         first call to routine COSQ1F or COSQ1B for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to COSQ1F and COSQ1B with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least N.


 Output Arguments
 
  R      Real output array.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(N-1)*INC).
 
         The output values of R are written over the input values.
         For J=0,...,N-1
 
          R(J*INC) =
 
              N-1
              SUM  R(N1*INC)*COS(J*(2*N1+1)*PI/(2*N))
              N1=0
 
 WSAVE   Contains values initialized by subroutine COSQ1I that
         must not be destroyed between calls to routine COSQ1F
         or COSQ1B.


 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSQ1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COSQ1F computes the one-dimensional Fourier 
 transform of a sequence which is a cosine series with odd wave 
 numbers.  This transform is referred to as the forward transform 
 or Fourier analysis, transforming the sequence from physical to 
 spectral space.
 
 This transform is normalized since a call to COSQ1F followed
 by a call to COSQ1B (or vice-versa) reproduces the original
 array subject to algorithmic constrain, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array with dimension LENSAV.  WSAVE's contents 
         must be initialized with a call to subroutine COSQ1I before 
         the first call to routine COSQ1F or COSQ1B for a given 
         transform length N.  WSAVE's contents may be re-used for 
         subsequent calls to COSQ1F and COSQ1B with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(N-1)*INC).
 
         The output values of R are written over the input values.
         For J=0,...,N-1
 
          R(J*INC) = 

              R(0)/N 
 
              N-1
            + SUM  2.*R(N1*INC)*COS((2*J+1)*N1*PI/(2*N))/N
              N1=1
 
 
 WSAVE   Contains values initialized by subroutine COSQ1I that
         must not be destroyed between calls to routine COSQ1F
         or COSQ1B.


 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSQMI (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine COSQMI initializes array WSAVE for use
 in its companion routines COSQMF and COSQMB.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines COSQMB or COSQMF.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COSQMB computes the one-dimensional Fourier 
 transform of multiple sequences, each of which is a cosine series 
 with odd wave numbers.  This transform is referred to as the 
 backward transform or Fourier synthesis, transforming the sequences 
 from spectral to physical space.
 
 This transform is normalized since a call to COSQMB followed
 by a call to COSQMF (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array with dimension LENSAV.  WSAVE's contents 
         must be initialized with a call to subroutine COSQMI before 
         the first call to routine COSQMF or COSQMB for a given 
         transform length N.  WSAVE's contents may be re-used for 
         subsequent calls to COSQMF and COSQMB with the same N.


 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*N.


 Output Arguments
 
  R      Real output array.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=0,...,N-1
 
          R(I*JUMP+J*INC) =
 
              N-1
              SUM  R(I*JUMP+N1*INC)*COS(J*(2*N1+1)*PI/(2*N))
              N1=0
 
 WSAVE   Contains values initialized by subroutine COSQMI that
         must not be destroyed between calls to routine COSQMF
         or COSQMB.

 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are 
              consistent if equality 
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N 
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables 
              INC, JUMP, N and LOT must be consistent, otherwise at 
              least one array element mistakenly is transformed more 
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE COSQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)

DESCRIPTION

 FFTPACK 5.1 routine COSQMF computes the one-dimensional Fourier 
 transform of multiple sequences within a real array, where each 
 of the sequences is a cosine series with odd wave numbers.  This 
 transform is referred to as the forward transform or Fourier 
 synthesis, transforming the sequences from spectral to physical 
 space.
 
 This transform is normalized since a call to COSQMF followed
 by a call to COSQMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.

 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array o length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine COSQMI before the 
         first call to routine COSQMF or COSQMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent
         calls to COSQMF and COSQMB with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real array of dimension LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(0:(LOT-1)*JUMP+(N-1)*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=0,...,N-1
 
          R(I*JUMP+J*INC) = 

              R(I*JUMP)/N 
 
              N-1
            + SUM  2.*R(I*JUMP+*N1*INC)*COS((2*J+1)*N1*PI/(2*N))/N
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are 
              consistent if equality 
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N 
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables 
              INC, JUMP, N and LOT must be consistent, otherwise at 
              least one array element mistakenly is transformed more 
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINQ1I (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine SINQ1I initializes array WSAVE for use
 in its companion routines SINQ1F and SINQ1B.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines SINQ1B or SINQ1F.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINQ1B (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINQ1B computes the one-dimensional Fourier 
 transform of a sequence which is a sine series with odd wave 
 numbers.  This transform is referred to as the backward transform 
 or Fourier synthesis, transforming the sequence from spectral to 
 physical space.
 
 This transform is normalized since a call to SINQ1B followed
 by a call to SINQ1F (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINQ1I before the 
         first call to routine SINQ1F or SINQ1B for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINQ1F and SINQ1B with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at least N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(INC:N*INC).
 
         The output values of R are written over the input values.
         For J=1,...,N
 
          R(J*INC) =
 
               N 
              SUM  R(N1*INC)*SIN(J*(2*N1-1)*PI/(2*N))
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINQ1F (N, INC, R, LENR, WSAVE, LENSAV, WORK, LENWRK, IER)

 INTEGER    N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINQ1F computes the one-dimensional Fourier 
 transform of a sequence which is a sine series of odd wave numbers.  
 This transform is referred to as the forward transform or Fourier 
 analysis, transforming the sequence from physical to spectral space.
 
 This transform is normalized since a call to SINQ1F followed
 by a call to SINQ1B (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 N       Integer length of the sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the sequence.
 
 R       Real array of length LENR containing the sequence to be 
         transformed.
 
 LENR    Integer dimension of R array.  LENR must be at least
         INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINQ1I before the 
         first call to routine SINQ1F or SINQ1B for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINQ1F and SINQ1B with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at least N.


 Output Arguments
 
 R       Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(INC:N*INC).
 
         The output values of R are written over the input values.
         For J=1,...,N
 
          R(J*INC) = 

              N-1
            + SUM  (2.*R(N1*INC)*SIN(((2*J-1)*N1*PI/(2*N)))/N
              N1=1

            + ((-1)**(J+1))*R(N*INC)/N
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINQMI (N, WSAVE, LENSAV, IER)
 INTEGER    N, LENSAV, IER
 REAL       WSAVE(LENSAV)
 

DESCRIPTION

 FFTPACK 5.1 subroutine SINQMI initializes array WSAVE for use
 in its companion routines SINQMF and SINQMB.  The prime factor-
 ization of N together with a tabulation of the trigonometric
 functions are computed and stored in array WSAVE.  Separate
 WSAVE arrays are required for different values of N.
 
 Input Arguments
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 Output Arguments
 
 WSAVE   Real work array with dimension LENSAV, containing the
         prime factors of N and also containing certain trigonometric 
         values which will be used in routines SINQMB or SINQMF.


 IER     Integer error return
         =  0 successful exit
         =  2 input parameter LENSAV not big enough
         = 20 input error returned by lower level routine

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINQMB (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINQMB computes the one-dimensional Fourier 
 transform of multiple sequences within a real array, where each 
 of the sequences is a sine series with odd wave numbers.  This 
 transform is referred to as the backward transform or Fourier 
 synthesis, transforming the sequences from spectral to physical 
 space.
 
 This transform is normalized since a call to SINQMB followed
 by a call to SINQMF (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINQMI before the 
         first call to routine SINQMF or SINQMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINQMF and SINQMB with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*N.


 Output Arguments
 
  R      Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(INC:(LOT-1)*JUMP+N*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=1,...,N
 
          R(I*JUMP+J*INC) =
 
              N 
              SUM  R(I*JUMP+N1*INC)*SIN(J*(2*N1-1)*PI/(2*N))
              N1=1
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are
              consistent if equality
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables
              INC, JUMP, N and LOT must be consistent ... otherwise at
              least one array element mistakenly is transformed more
              than once.

C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C     *                                                               *
C     *                  copyright (c) 2011 by UCAR                   *
C     *                                                               *
C     *       University Corporation for Atmospheric Research         *
C     *                                                               *
C     *                      all rights reserved                      *
C     *                                                               *
C     *                     FFTPACK  version 5.1                      *
C     *                                                               *
C     *                 A Fortran Package of Fast Fourier             *
C     *                                                               *
C     *                Subroutines and Example Programs               *
C     *                                                               *
C     *                             by                                *
C     *                                                               *
C     *               Paul Swarztrauber and Dick Valent               *
C     *                                                               *
C     *                             of                                *
C     *                                                               *
C     *         the National Center for Atmospheric Research          *
C     *                                                               *
C     *                Boulder, Colorado  (80307)  U.S.A.             *
C     *                                                               *
C     *                   which is sponsored by                       *
C     *                                                               *
C     *              the National Science Foundation                  *
C     *                                                               *
C     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *

SYNOPSIS

 SUBROUTINE SINQMF (LOT, JUMP, N, INC, R, LENR, WSAVE, LENSAV, 
1                   WORK, LENWRK, IER)

 INTEGER    LOT, JUMP, N, INC, LENR, LENSAV, LENWRK, IER
 REAL       R(LENR), WSAVE(LENSAV), WORK(LENWRK)
 

DESCRIPTION

 FFTPACK 5.1 routine SINQMF computes the one-dimensional Fourier 
 transform of multiple sequences within a real array, where each 
 sequence is a sine series with odd wave numbers.  This transform 
 is referred to as the forward transform or Fourier synthesis, 
 transforming the sequences from spectral to physical space.
 
 This transform is normalized since a call to SINQMF followed
 by a call to SINQMB (or vice-versa) reproduces the original
 array subject to algorithmic constraints, roundoff error, etc.
 
 Input Arguments
 
 LOT     Integer number of sequences to be transformed within
         array R.
 
 JUMP    Integer increment between the locations, in array R,
         of the first elements of two consecutive sequences
         to be transformed.
 
 N       Integer length of each sequence to be transformed.  The 
         transform is most efficient when N is a product of 
         small primes.
 
 INC     Integer increment between the locations, in array R,
         of two consecutive elements within the same sequence.
 
 R       Real array containing LOT sequences, each having length N.
         R can have any number of dimensions, but the total number 
         of locations must be at least LENR.
 
 LENR    Integer dimension of R array.  LENR must be at least
         (LOT-1)*JUMP + INC*(N-1)+ 1.


 WSAVE   Real work array of length LENSAV.  WSAVE's contents must 
         be initialized with a call to subroutine SINQMI before the 
         first call to routine SINQMF or SINQMB for a given transform
         length N.  WSAVE's contents may be re-used for subsequent 
         calls to SINQMF and SINQMB with the same N.
 
 LENSAV  Integer dimension of WSAVE array.  LENSAV must be at least 
         2*N + INT(LOG (REAL(N))/LOG(2.)) +4.


 WORK    Real work array of dimension at least LENWRK.
 
 LENWRK  Integer dimension of WORK array.  LENWRK must be at 
         least LOT*N.


 Output Arguments
 
 R       Real output array R.  For purposes of exposition, 
         assume R's range of indices is given by 
         R(INC:(LOT-1)*JUMP+N*INC).
 
         The output values of R are written over the input values.
         For I=0,...,LOT-1 and J=1,...,N
 
          R(I*JUMP+J*INC) = 

              N-1
            + SUM  (2.*R(I*JUMP+*N1*INC)*SIN(((2*J-1)*N1*PI/(2*N)))/N
              N1=1

            + ((-1)**(J+1))*R(I*JUMP+N*INC)/N
 
 IER     Integer error return
         =  0 successful exit
         =  1 input parameter LENR   not big enough
         =  2 input parameter LENSAV not big enough
         =  3 input parameter LENWRK not big enough
         =  4 input parameters INC,JUMP,N,LOT are not consistent.
         = 20 input error returned by lower level routine

              The parameters integers INC, JUMP, N and LOT are
              consistent if equality
              I1*INC + J1*JUMP = I2*INC + J2*JUMP for I1,I2 < N
              and J1,J2 < LOT implies I1=I2 and J1=J2.

              For multiple FFTs to execute correctly, input variables
              INC, JUMP, N and LOT must be consistent ... otherwise at
              least one array element mistakenly is transformed more
              than once.