FFTPACK5 Documentation

 

 

FFTPACK 5.1 -- a FORTRAN library of fast Fourier transforms

Complex Transform Routines
  • CFFTMI multiple complex initialization
  • CFFTMB multiple complex backward
  • CFFTMF multiple complex forward
Real Transform Routines
  • 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

CFFT1I - initialization routine for CFFT1B and CFFT1F

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


CFFT1B - complex backward fast Fourier transform

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, 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


CFFT1F - complex forward fast Fourier transform

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, 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


CFFT2I - initialization routine for CFFT2B, CFFT2F

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


CFFT2B - complex, two-dimensional backward fast Fourier transform

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, 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


CFFT2F - complex, two-dimensional forward fast Fourier transform

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, 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


CFFTMI - initialization routine for CFFTMB and CFFTMF

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


CFFTMB - complex, multiple backward fast Fourier transform

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, 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.


CFFTMF - complex, multiple forward fast Fourier transform

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, 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.


RFFT1I - initialization routine for RFFT1B and RFFT1F

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


RFFT1B - real backward fast Fourier transform

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


RFFT1F - real backward fast Fourier transform

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


RFFT2I - initialization routine for RFFT2B and RFFT2F

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


RFFT2B - complex to real, two-dimensional
backward fast Fourier transform

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


RFFT2F - real to complex, two-dimensional
forward fast Fourier transform

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


RFFTMI - initialization routine for RFFTMB and RFFTMF

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


RFFTMB - real, multiple backward fast Fourier transform

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, 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.


RFFTMF - real, multiple forward fast Fourier transform

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, 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.


COST1I - initialization routine for COST1B and COST1F

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


COST1B - real backward cosine fast Fourier transform

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


COST1F - real backward cosine fast Fourier transform

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


COSTMI - initialization routine for COSTMB and COSTMF

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


COSTMB - real, multiple backward cosine fast Fourier transform

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, 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.


COSTMF - real, multiple forward cosine fast Fourier transform

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, 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.


SINT1I - initialization routine for SINT1B and SINT1F

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


SINT1B - real backward sine fast Fourier transform

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


SINT1F - real forward sine fast Fourier transform

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


SINTMI - initialization routine for SINTMB and SINTMF

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



SINTMB - real, multiple backward sine fast Fourier transform

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, 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.


SINTMF - real, multiple forward sine fast Fourier transform

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, 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.


COSQ1I - initialization routine for COSQ1B and COSQ1F

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

COSQ1B - real, backward quarter-cosine fast Fourier transform

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


COSQ1F - real, forward quarter-cosine fast Fourier transform

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

COSQMI - initialization routine for COSQMB and COSQMF

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


COSQMB - real, multiple backward quarter-cosine fast Fourier transform

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, 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.


COSQMF - real, multiple forward quarter-cosine fast Fourier transform

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, 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.


SINQ1I - initialization routine for SINQ1B and SINQ1F

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


SINQ1B - real backward quarter-sine fast Fourier transform

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


SINQ1F - real forward quarter-sine fast Fourier transform

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


SINQMI - initialization routine for SINQMB and SINQMF

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


SINQMB - real, multiple backward quarter-sine fast Fourier transform

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, 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.


SINQMF - real, multiple forward quarter-sine fast Fourier transform

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, 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.

FFTPACK5, Copyright (C) 2004-2011, Computational Information Systems Laboratory, University Corporation for Atmospheric Research