c
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c  .                                                             .
c  .                  copyright (c) 1998 by UCAR                 .
c  .                                                             .
c  .       University Corporation for Atmospheric Research       .
c  .                                                             .
c  .                      all rights reserved                    .
c  .                                                             .
c  .                                                             .
c  .                         SPHEREPACK                          .
c  .                                                             .
c  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
c
c
c ... file hrfft.f
c
c     this file contains a multiple fft package for spherepack3.0.
c     it includes code and documentation for performing fast fourier
c     transforms (see subroutines hrffti,hrfftf and hrfftb)
c
c **********************************************************************
c
c     subroutine hrffti(n,wsave)
c
c     subroutine hrffti initializes the array wsave which is used in
c     both hrfftf and hrfftb. the prime factorization of n together 
c     with a tabulation of the trigonometric functions are computed and
c     stored in wsave.
c
c     input parameter
c
c     n       the length of the sequence to be transformed.
c
c     output parameter
c
c     wsave   a work array which must be dimensioned at least 2*n+15.
c             the same work array can be used for both hrfftf and 
c             hrfftb as long as n remains unchanged. different wsave 
c             arrays are required for different values of n. the 
c             contents of wsave must not be changed between calls 
c             of hrfftf or hrfftb.
c
c **********************************************************************
c
c     subroutine hrfftf(m,n,r,mdimr,wsave,work)
c
c     subroutine hrfftf computes the fourier coefficients of m real
c     perodic sequences (fourier analysis); i.e. hrfftf computes the
c     real fft of m sequences each with length n. the transform is 
c     defined below at output parameter r.
c
c     input parameters
c
c     m       the number of sequences.
c
c     n       the length of all m sequences.  the method is most
c             efficient when n is a product of small primes. n may
c             change as long as different work arrays are provided
c
c     r       r(m,n) is a two dimensional real array that contains m
c             sequences each with length n.
c
c     mdimr   the first dimension of the r array as it appears
c             in the program that calls hrfftf. mdimr must be
c             greater than or equal to m.
c
c
c     wsave   a work array with at least least 2*n+15 locations
c             in the program that calls hrfftf. the wsave array must be
c             initialized by calling subroutine hrffti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c             the same wsave array can be used by hrfftf and hrfftb.
c
c     work    a real work array with m*n locations.
c
c
c     output parameters
c
c     r      for all j=1,...,m
c  
c             r(j,1) = the sum from i=1 to i=n of r(j,i)
c
c             if n is even set l =n/2   , if n is odd set l = (n+1)/2
c
c               then for k = 2,...,l
c
c                  r(j,2*k-2) = the sum from i = 1 to i = n of
c
c                       r(j,i)*cos((k-1)*(i-1)*2*pi/n)
c
c                  r(j,2*k-1) = the sum from i = 1 to i = n of
c
c                      -r(j,i)*sin((k-1)*(i-1)*2*pi/n)
c
c             if n is even
c
c                  r(j,n) = the sum from i = 1 to i = n of
c
c                       (-1)**(i-1)*r(j,i)
c
c      *****  note
c                  this transform is unnormalized since a call of hrfftf
c                  followed by a call of hrfftb will multiply the input
c                  sequence by n.
c
c     wsave   contains results which must not be destroyed between
c             calls of hrfftf or hrfftb.
c
c     work    a real work array with m*n locations that does
c             not have to be saved.
c
c **********************************************************************
c
c     subroutine hrfftb(m,n,r,mdimr,wsave,work)
c
c     subroutine hrfftb computes the real perodic sequence of m
c     sequences from their fourier coefficients (fourier synthesis). 
c     the transform is defined below at output parameter r.
c
c     input parameters
c
c     m       the number of sequences.
c
c     n       the length of all m sequences.  the method is most
c             efficient when n is a product of small primes. n may
c             change as long as different work arrays are provided
c
c     r       r(m,n) is a two dimensional real array that contains
c             the fourier coefficients of m sequences each with 
c             length n.
c
c     mdimr   the first dimension of the r array as it appears
c             in the program that calls hrfftb. mdimr must be
c             greater than or equal to m.
c
c     wsave   a work array which must be dimensioned at least 2*n+15.
c             in the program that calls hrfftb. the wsave array must be
c             initialized by calling subroutine hrffti(n,wsave) and a
c             different wsave array must be used for each different
c             value of n. this initialization does not have to be
c             repeated so long as n remains unchanged thus subsequent
c             transforms can be obtained faster than the first.
c             the same wsave array can be used by hrfftf and hrfftb.
c
c     work    a real work array with m*n locations.
c
c
c     output parameters
c
c     r      for all j=1,...,m
c  
c             for n even and for i = 1,...,n
c
c                  r(j,i) = r(j,1)+(-1)**(i-1)*r(j,n)
c
c                       plus the sum from k=2 to k=n/2 of
c
c                        2.*r(j,2*k-2)*cos((k-1)*(i-1)*2*pi/n)
c
c                       -2.*r(j,2*k-1)*sin((k-1)*(i-1)*2*pi/n)
c
c             for n odd and for i = 1,...,n
c
c                  r(j,i) = r(j,1) plus the sum from k=2 to k=(n+1)/2 of
c
c                       2.*r(j,2*k-2)*cos((k-1)*(i-1)*2*pi/n)
c
c                      -2.*r(j,2*k-1)*sin((k-1)*(i-1)*2*pi/n)
c
c      *****  note
c                  this transform is unnormalized since a call of hrfftf
c                  followed by a call of hrfftb will multiply the input
c                  sequence by n.
c
c     wsave   contains results which must not be destroyed between
c             calls of hrfftb or hrfftf.
c
c     work    a real work array with m*n locations that does not
c             have to be saved
c
c **********************************************************************
c
c
c