FUNCTION pimach(dum)
!***BEGIN PROLOGUE  PIMACH
!
!  This subprogram supplies the value of the constant PI correct to
!  machine precision where
!
!  PI=3.1415926535897932384626433832795028841971693993751058209749446
!***ROUTINES CALLED  (NONE)
!***END PROLOGUE  PIMACH
!
!***FIRST EXECUTABLE STATEMENT  PIMACH
  pimach = 3.14159265358979
  RETURN
  END FUNCTION pimach



SUBROUTINE vcost(m,n,x,xt,mdimx,wsave) 8,1

!***BEGIN PROLOGUE  VCOST
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  900509   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE
!          SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NIST)
!***PURPOSE  Cosine transform of one or more real, even sequences.
!***DESCRIPTION
!
!  Subroutine VCOST computes the discrete Fourier cosine transform
!  of M even sequences X(J,I), J=1,...,M.  The transform is defined
!  below at output parameter X.
!
!  The array WSAVE which is used by subroutine VCOST must be
!  initialized by calling subroutine VCOSTI(N,WSAVE).
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequence to be transformed.  N must be
!       greater than 1.  The method is most efficient when N-1 is
!       is a product of small primes.
!
!  X       an array of size at least X(MDIMX,N) which contains the
!       the sequences to be transformed.  The sequences are stored
!       in the ROWS of X.  Thus, the Jth sequence is stored in
!       X(J,I), I=1,..,N.
!
!  XT      a work array of size at least XT(MDIMX,N-1).
!
!  MDIMX   the first dimension of the array X exactly as it appears in
!       the calling program.
!
!  WSAVE   a work array which must be dimensioned at least 3*N+15
!       in the program that calls VCOST.  The WSAVE array must be
!       initialized by calling subroutine VCOSTI(N,WSAVE), and a
!       different WSAVE array must be used for each different
!       value of N.  This initialization does not have to be
!       repeated so long as N remains unchanged.  Thus subsequent
!       transforms can be obtained faster than the first.
!
!  Output Parameters
!
!  X       For I=1,...,N and J=1,...,M
!
!          X(J,I) = ( X(J,1)+(-1)**(I-1)*X(J,N)
!
!            + the sum from K=2 to K=N-1
!
!              2*X(J,K)*COS((K-1)*(I-1)*PI/(N-1)) )/SQRT(2*(N-1))
!
!  WSAVE   contains initialization calculations which must not be
!       destroyed between calls of VCOST.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VCOST followed immediately by another call
!        of VCOST will return the original sequences X.  Thus,
!        VCOST is the correctly normalized inverse of itself.
!
!  -----------------------------------------------------------------
!
!  VCOST is a straightforward extension of the subprogram COST to
!  handle M simultaneous sequences.  The scaling of the sequences
!  computed by VCOST is different than that of COST.  COST was
!  originally developed by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!            Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!            pp. 51-83.
!***ROUTINES CALLED  VRFFTF
!***END PROLOGUE  VCOST
  DIMENSION       x(mdimx,*), xt(mdimx,*), wsave(*)
!***FIRST EXECUTABLE STATEMENT  VCOST
  IF (m <= 0)  GO TO 900
  IF (n <= 1)  GO TO 900
  IF (n > 3)  GO TO 400
  IF (n == 3)  GO TO 300
!
!  CASE  N = 2
!
  scale = SQRT(0.50E0)
  DO j=1,m
    x1h = scale*(x(j,1)+x(j,2))
    x(j,2) = scale*(x(j,1)-x(j,2))
    x(j,1) = x1h
  END DO
  GO TO 900
!
!  CASE  N = 3
!
  300 CONTINUE
  scale = 0.50E0
  DO j=1,m
    x1p3 = x(j,1)+x(j,3)
    tx2 = x(j,2)+x(j,2)
    x(j,2) = scale*(x(j,1)-x(j,3))
    x(j,1) = scale*(x1p3+tx2)
    x(j,3) = scale*(x1p3-tx2)
  END DO
  GO TO 900
!
!  CASE  N .GT. 3
!
!  ... PREPROCESSING
!
  400 CONTINUE
  nm1 = n-1
  np1 = n+1
  ns2 = n/2
  DO j=1,m
    xt(j,1) = x(j,1)-x(j,n)
    x(j,1) = x(j,1)+x(j,n)
  END DO
  DO k=2,ns2
    kc = np1-k
    DO j=1,m
      t1 = x(j,k)+x(j,kc)
      t2 = x(j,k)-x(j,kc)
      xt(j,1) = xt(j,1)+wsave(kc)*t2
      t2 = wsave(k)*t2
      x(j,k) = t1-t2
      x(j,kc) = t1+t2
    END DO
  END DO
  modn = MOD(n,2)
  IF (modn /= 0) THEN
    DO j=1,m
      x(j,ns2+1) = x(j,ns2+1)+x(j,ns2+1)
    END DO
  END IF
  DO j=1,m
    x(j,n) = xt(j,1)
  END DO
!
!  ... REAL PERIODIC TRANSFORM
!
  CALL vrfftf (m,nm1,x,xt,mdimx,wsave(np1))
!
!  ... POSTPROCESSING
!
  factor = 1.0/SQRT(REAL(nm1))
  DO j=1,m
    xt(j,1) = x(j,2)
    x(j,2) = factor*x(j,n)
  END DO
  DO i=4,n,2
    DO j=1,m
      xi = x(j,i)
      x(j,i) = x(j,i-2)-x(j,i-1)
      x(j,i-1) = xt(j,1)
      xt(j,1) = xi
    END DO
  END DO
  IF (modn /= 0) THEN
    DO j=1,m
      x(j,n) = xt(j,1)
    END DO
  END IF
!
!  ... NORMALIZATION
!
  scale = SQRT(0.5)
  DO i=1,n
    DO j=1,m
      x(j,i) = scale*x(j,i)
    END DO
  END DO
!
!  EXIT
!
  900 CONTINUE
  RETURN
END SUBROUTINE vcost


SUBROUTINE vcosti(n,wsave) 8,1
!***BEGIN PROLOGUE  VCOSTI
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  900509   (YYMMDD)
!***CATEGORY NO.  J1A3
!***KEYWORDS  FAST FOURIER TRANSFORM, COSINE TRANSFORM, MULTIPLE
!          SEQUENCES
!***AUTHOR  BOISVERT, R. F. (NIST)
!***PURPOSE  Initialize for VCOST.
!***DESCRIPTION
!
!  Subroutine VCOSTI initializes the array WSAVE which is used in
!  subroutine VCOST.  The prime factorization of N together with
!  a tabulation of the trigonometric functions are computed and
!  stored in WSAVE.
!
!  Input Parameter
!
!  N       the length of the sequence to be transformed.  The method
!       is most efficient when N-1 is a product of small primes.
!
!  Output Parameter
!
!  WSAVE   a work array which must be dimensioned at least 3*N+15.
!       Different WSAVE arrays are required for different values
!       of N.  The contents of WSAVE must not be changed between
!       calls of VCOST.
!
!  -----------------------------------------------------------------
!
!  VCOSTI is a straightforward extension of the subprogram COSTI to
!  handle M simultaneous sequences.  COSTI was originally developed
!  by P. N. Swarztrauber of NCAR.
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!            Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!            pp. 51-83.
!***ROUTINES CALLED  VRFFTI
!***END PROLOGUE  VCOSTI
  DIMENSION       wsave(*)
!***FIRST EXECUTABLE STATEMENT  VCOSTI
  pi = pimach(1.0)
  IF (n <= 3) RETURN
  nm1 = n-1
  np1 = n+1
  ns2 = n/2
  dt = pi/REAL(nm1)
  fk = 0.
  DO k=2,ns2
    fk = fk+1.
    wsave(k) = 2.*SIN(fk*dt)
  END DO
  fk = 0.
  DO k=2,ns2
    kc = np1-k
    fk = fk+1.
    wsave(kc) = 2.*COS(fk*dt)
  END DO
  CALL vrffti (nm1,wsave(n+1))
  RETURN
END SUBROUTINE vcosti


SUBROUTINE vradb2 (mp,ido,l1,cc,ch,mdimc,wa1) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION  cc(mdimc,ido,2,l1)    ,ch(mdimc,ido,l1,2),                 &
                  wa1(ido)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,k,1) = cc(m,1,1,k)+cc(m,ido,2,k)
      ch(m,1,k,2) = cc(m,1,1,k)-cc(m,ido,2,k)
    END DO
  END DO
  IF (ido-2 < 0) THEN
    GO TO   107
  ELSE IF (ido-2 == 0) THEN
    GO TO   105
  END IF
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,k,1) = cc(m,i-1,1,k)+cc(m,ic-1,2,k)
        ch(m,i,k,1) = cc(m,i,1,k)-cc(m,ic,2,k)
        ch(m,i-1,k,2) = wa1(i-2)*(cc(m,i-1,1,k)-cc(m,ic-1,2,k))         &
            -wa1(i-1)*(cc(m,i,1,k)+cc(m,ic,2,k))
        ch(m,i,k,2) = wa1(i-2)*(cc(m,i,1,k)+cc(m,ic,2,k))+wa1(i-1)      &
            *(cc(m,i-1,1,k)-cc(m,ic-1,2,k))
      END DO
    END DO
  END DO
  IF (MOD(ido,2) == 1) RETURN
  105 DO k=1,l1
    DO m=1,mp
      ch(m,ido,k,1) = cc(m,ido,1,k)+cc(m,ido,1,k)
      ch(m,ido,k,2) = -(cc(m,1,2,k)+cc(m,1,2,k))
    END DO
  END DO
  107 RETURN
END SUBROUTINE vradb2


SUBROUTINE vradb3 (mp,ido,l1,cc,ch,mdimc,wa1,wa2) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION  cc(mdimc,ido,3,l1)    ,ch(mdimc,ido,l1,3),                 &
                  wa1(ido)   ,wa2(ido)
  arg=2.*pimach(1.0)/3.
  taur=COS(arg)
  taui=SIN(arg)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,k,1) = cc(m,1,1,k)+2.*cc(m,ido,2,k)
      ch(m,1,k,2) = cc(m,1,1,k)+(2.*taur)*cc(m,ido,2,k)                 &
          -(2.*taui)*cc(m,1,3,k)
      ch(m,1,k,3) = cc(m,1,1,k)+(2.*taur)*cc(m,ido,2,k)                 &
          +2.*taui*cc(m,1,3,k)
    END DO
  END DO
  IF (ido == 1) RETURN
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,k,1) = cc(m,i-1,1,k)+(cc(m,i-1,3,k)+cc(m,ic-1,2,k))
        ch(m,i,k,1) = cc(m,i,1,k)+(cc(m,i,3,k)-cc(m,ic,2,k))
        ch(m,i-1,k,2) = wa1(i-2)*                                       &
            ((cc(m,i-1,1,k)+taur*(cc(m,i-1,3,k)+cc(m,ic-1,2,k)))-       &
            (taui*(cc(m,i,3,k)+cc(m,ic,2,k))))                          &
                     -wa1(i-1)*                                         &
            ((cc(m,i,1,k)+taur*(cc(m,i,3,k)-cc(m,ic,2,k)))+             &
            (taui*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))))
        ch(m,i,k,2) = wa1(i-2)*                                         &
            ((cc(m,i,1,k)+taur*(cc(m,i,3,k)-cc(m,ic,2,k)))+             &
            (taui*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))))                      &
                    +wa1(i-1)*                                          &
            ((cc(m,i-1,1,k)+taur*(cc(m,i-1,3,k)+cc(m,ic-1,2,k)))-       &
            (taui*(cc(m,i,3,k)+cc(m,ic,2,k))))
        ch(m,i-1,k,3) = wa2(i-2)*                                       &
            ((cc(m,i-1,1,k)+taur*(cc(m,i-1,3,k)+cc(m,ic-1,2,k)))+       &
            (taui*(cc(m,i,3,k)+cc(m,ic,2,k))))                          &
                      -wa2(i-1)*                                        &
            ((cc(m,i,1,k)+taur*(cc(m,i,3,k)-cc(m,ic,2,k)))-             &
            (taui*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))))
        ch(m,i,k,3) = wa2(i-2)*                                         &
            ((cc(m,i,1,k)+taur*(cc(m,i,3,k)-cc(m,ic,2,k)))-             &
            (taui*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))))                      &
                   +wa2(i-1)*                                           &
            ((cc(m,i-1,1,k)+taur*(cc(m,i-1,3,k)+cc(m,ic-1,2,k)))+       &
            (taui*(cc(m,i,3,k)+cc(m,ic,2,k))))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE vradb3


SUBROUTINE vradb4 (mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION  cc(mdimc,ido,4,l1)  ,ch(mdimc,ido,l1,4)    ,               &
                  wa1(ido)  ,wa2(ido)  ,wa3(ido)
  sqrt2=SQRT(2.)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,k,3) = (cc(m,1,1,k)+cc(m,ido,4,k))                         &
          -(cc(m,ido,2,k)+cc(m,ido,2,k))
      ch(m,1,k,1) = (cc(m,1,1,k)+cc(m,ido,4,k))                         &
          +(cc(m,ido,2,k)+cc(m,ido,2,k))
      ch(m,1,k,4) = (cc(m,1,1,k)-cc(m,ido,4,k))                         &
          +(cc(m,1,3,k)+cc(m,1,3,k))
      ch(m,1,k,2) = (cc(m,1,1,k)-cc(m,ido,4,k))                         &
          -(cc(m,1,3,k)+cc(m,1,3,k))
    END DO
  END DO
  IF (ido-2 < 0) THEN
    GO TO   107
  ELSE IF (ido-2 == 0) THEN
    GO TO   105
  END IF
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,k,1) = (cc(m,i-1,1,k)+cc(m,ic-1,4,k))                  &
            +(cc(m,i-1,3,k)+cc(m,ic-1,2,k))
        ch(m,i,k,1) = (cc(m,i,1,k)-cc(m,ic,4,k))                        &
            +(cc(m,i,3,k)-cc(m,ic,2,k))
        ch(m,i-1,k,2)=wa1(i-2)*((cc(m,i-1,1,k)-cc(m,ic-1,4,k))          &
            -(cc(m,i,3,k)+cc(m,ic,2,k)))-wa1(i-1)                       &
            *((cc(m,i,1,k)+cc(m,ic,4,k))+(cc(m,i-1,3,k)-cc(m,ic-1,2,k)))
        ch(m,i,k,2)=wa1(i-2)*((cc(m,i,1,k)+cc(m,ic,4,k))                &
            +(cc(m,i-1,3,k)-cc(m,ic-1,2,k)))+wa1(i-1)                   &
            *((cc(m,i-1,1,k)-cc(m,ic-1,4,k))-(cc(m,i,3,k)+cc(m,ic,2,k)))
        ch(m,i-1,k,3)=wa2(i-2)*((cc(m,i-1,1,k)+cc(m,ic-1,4,k))          &
            -(cc(m,i-1,3,k)+cc(m,ic-1,2,k)))-wa2(i-1)                   &
            *((cc(m,i,1,k)-cc(m,ic,4,k))-(cc(m,i,3,k)-cc(m,ic,2,k)))
        ch(m,i,k,3)=wa2(i-2)*((cc(m,i,1,k)-cc(m,ic,4,k))                &
            -(cc(m,i,3,k)-cc(m,ic,2,k)))+wa2(i-1)                       &
            *((cc(m,i-1,1,k)+cc(m,ic-1,4,k))-(cc(m,i-1,3,k)             &
            +cc(m,ic-1,2,k)))
        ch(m,i-1,k,4)=wa3(i-2)*((cc(m,i-1,1,k)-cc(m,ic-1,4,k))          &
            +(cc(m,i,3,k)+cc(m,ic,2,k)))-wa3(i-1)                       &
            *((cc(m,i,1,k)+cc(m,ic,4,k))-(cc(m,i-1,3,k)-cc(m,ic-1,2,k)))
        ch(m,i,k,4)=wa3(i-2)*((cc(m,i,1,k)+cc(m,ic,4,k))                &
            -(cc(m,i-1,3,k)-cc(m,ic-1,2,k)))+wa3(i-1)                   &
            *((cc(m,i-1,1,k)-cc(m,ic-1,4,k))+(cc(m,i,3,k)+cc(m,ic,2,k)))
      END DO
    END DO
  END DO
  IF (MOD(ido,2) == 1) RETURN
  105 CONTINUE
  DO k=1,l1
    DO m=1,mp
      ch(m,ido,k,1) = (cc(m,ido,1,k)+cc(m,ido,3,k))                     &
          +(cc(m,ido,1,k)+cc(m,ido,3,k))
      ch(m,ido,k,2) = sqrt2*((cc(m,ido,1,k)-cc(m,ido,3,k))              &
          -(cc(m,1,2,k)+cc(m,1,4,k)))
      ch(m,ido,k,3) = (cc(m,1,4,k)-cc(m,1,2,k))                         &
          +(cc(m,1,4,k)-cc(m,1,2,k))
      ch(m,ido,k,4) = -sqrt2*((cc(m,ido,1,k)-cc(m,ido,3,k))             &
          +(cc(m,1,2,k)+cc(m,1,4,k)))
    END DO
  END DO
  107 RETURN
END SUBROUTINE vradb4


SUBROUTINE vradb5 (mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3,wa4) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION  cc(mdimc,ido,5,l1)    ,ch(mdimc,ido,l1,5),                 &
               wa1(ido)     ,wa2(ido)     ,wa3(ido)     ,wa4(ido)
  arg=2.*pimach(1.0)/5.
  tr11=COS(arg)
  ti11=SIN(arg)
  tr12=COS(2.*arg)
  ti12=SIN(2.*arg)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,k,1) = cc(m,1,1,k)+2.*cc(m,ido,2,k)+2.*cc(m,ido,4,k)
      ch(m,1,k,2) = (cc(m,1,1,k)+tr11*2.*cc(m,ido,2,k)                  &
          +tr12*2.*cc(m,ido,4,k))-(ti11*2.*cc(m,1,3,k)                  &
          +ti12*2.*cc(m,1,5,k))
      ch(m,1,k,3) = (cc(m,1,1,k)+tr12*2.*cc(m,ido,2,k)                  &
          +tr11*2.*cc(m,ido,4,k))-(ti12*2.*cc(m,1,3,k)                  &
          -ti11*2.*cc(m,1,5,k))
      ch(m,1,k,4) = (cc(m,1,1,k)+tr12*2.*cc(m,ido,2,k)                  &
          +tr11*2.*cc(m,ido,4,k))+(ti12*2.*cc(m,1,3,k)                  &
          -ti11*2.*cc(m,1,5,k))
      ch(m,1,k,5) = (cc(m,1,1,k)+tr11*2.*cc(m,ido,2,k)                  &
          +tr12*2.*cc(m,ido,4,k))+(ti11*2.*cc(m,1,3,k)                  &
          +ti12*2.*cc(m,1,5,k))
    END DO
  END DO
  IF (ido == 1) RETURN
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,k,1) = cc(m,i-1,1,k)+(cc(m,i-1,3,k)+cc(m,ic-1,2,k))    &
            +(cc(m,i-1,5,k)+cc(m,ic-1,4,k))
        ch(m,i,k,1) = cc(m,i,1,k)+(cc(m,i,3,k)-cc(m,ic,2,k))            &
            +(cc(m,i,5,k)-cc(m,ic,4,k))
        ch(m,i-1,k,2) = wa1(i-2)*((cc(m,i-1,1,k)+tr11*                  &
            (cc(m,i-1,3,k)+cc(m,ic-1,2,k))+tr12                         &
            *(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))-(ti11*(cc(m,i,3,k)         &
            +cc(m,ic,2,k))+ti12*(cc(m,i,5,k)+cc(m,ic,4,k))))            &
            -wa1(i-1)*((cc(m,i,1,k)+tr11*(cc(m,i,3,k)-cc(m,ic,2,k))     &
            +tr12*(cc(m,i,5,k)-cc(m,ic,4,k)))+(ti11*(cc(m,i-1,3,k)      &
            -cc(m,ic-1,2,k))+ti12*(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))
        ch(m,i,k,2) = wa1(i-2)*((cc(m,i,1,k)+tr11*(cc(m,i,3,k)          &
            -cc(m,ic,2,k))+tr12*(cc(m,i,5,k)-cc(m,ic,4,k)))             &
            +(ti11*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))+ti12                  &
            *(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))+wa1(i-1)                  &
            *((cc(m,i-1,1,k)+tr11*(cc(m,i-1,3,k)                        &
            +cc(m,ic-1,2,k))+tr12*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))       &
            -(ti11*(cc(m,i,3,k)+cc(m,ic,2,k))+ti12                      &
            *(cc(m,i,5,k)+cc(m,ic,4,k))))
        ch(m,i-1,k,3) = wa2(i-2)                                        &
            *((cc(m,i-1,1,k)+tr12*(cc(m,i-1,3,k)+cc(m,ic-1,2,k))        &
            +tr11*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))-(ti12*(cc(m,i,3,k)    &
            +cc(m,ic,2,k))-ti11*(cc(m,i,5,k)+cc(m,ic,4,k))))            &
            -wa2(i-1)                                                   &
            *((cc(m,i,1,k)+tr12*(cc(m,i,3,k)-                           &
            cc(m,ic,2,k))+tr11*(cc(m,i,5,k)-cc(m,ic,4,k)))              &
            +(ti12*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))-ti11                  &
            *(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))
        ch(m,i,k,3) = wa2(i-2)                                          &
            *((cc(m,i,1,k)+tr12*(cc(m,i,3,k)-                           &
            cc(m,ic,2,k))+tr11*(cc(m,i,5,k)-cc(m,ic,4,k)))              &
            +(ti12*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))-ti11                  &
            *(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))                           &
            +wa2(i-1)                                                   &
            *((cc(m,i-1,1,k)+tr12*(cc(m,i-1,3,k)+cc(m,ic-1,2,k))        &
            +tr11*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))-(ti12*(cc(m,i,3,k)    &
            +cc(m,ic,2,k))-ti11*(cc(m,i,5,k)+cc(m,ic,4,k))))
        ch(m,i-1,k,4) = wa3(i-2)                                        &
            *((cc(m,i-1,1,k)+tr12*(cc(m,i-1,3,k)+cc(m,ic-1,2,k))        &
            +tr11*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+(ti12*(cc(m,i,3,k)    &
            +cc(m,ic,2,k))-ti11*(cc(m,i,5,k)+cc(m,ic,4,k))))            &
            -wa3(i-1)                                                   &
            *((cc(m,i,1,k)+tr12*(cc(m,i,3,k)-                           &
            cc(m,ic,2,k))+tr11*(cc(m,i,5,k)-cc(m,ic,4,k)))              &
            -(ti12*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))-ti11                  &
            *(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))
        ch(m,i,k,4) = wa3(i-2)                                          &
            *((cc(m,i,1,k)+tr12*(cc(m,i,3,k)-                           &
            cc(m,ic,2,k))+tr11*(cc(m,i,5,k)-cc(m,ic,4,k)))              &
            -(ti12*(cc(m,i-1,3,k)-cc(m,ic-1,2,k))-ti11                  &
            *(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))                           &
            +wa3(i-1)                                                   &
            *((cc(m,i-1,1,k)+tr12*(cc(m,i-1,3,k)+cc(m,ic-1,2,k))        &
            +tr11*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+(ti12*(cc(m,i,3,k)    &
            +cc(m,ic,2,k))-ti11*(cc(m,i,5,k)+cc(m,ic,4,k))))
        ch(m,i-1,k,5) = wa4(i-2)                                        &
            *((cc(m,i-1,1,k)+tr11*(cc(m,i-1,3,k)+cc(m,ic-1,2,k))        &
            +tr12*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+(ti11*(cc(m,i,3,k)    &
            +cc(m,ic,2,k))+ti12*(cc(m,i,5,k)+cc(m,ic,4,k))))            &
            -wa4(i-1)                                                   &
            *((cc(m,i,1,k)+tr11*(cc(m,i,3,k)-cc(m,ic,2,k))              &
            +tr12*(cc(m,i,5,k)-cc(m,ic,4,k)))-(ti11*(cc(m,i-1,3,k)      &
            -cc(m,ic-1,2,k))+ti12*(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))
        ch(m,i,k,5) = wa4(i-2)                                          &
            *((cc(m,i,1,k)+tr11*(cc(m,i,3,k)-cc(m,ic,2,k))              &
            +tr12*(cc(m,i,5,k)-cc(m,ic,4,k)))-(ti11*(cc(m,i-1,3,k)      &
            -cc(m,ic-1,2,k))+ti12*(cc(m,i-1,5,k)-cc(m,ic-1,4,k))))      &
            +wa4(i-1)                                                   &
            *((cc(m,i-1,1,k)+tr11*(cc(m,i-1,3,k)+cc(m,ic-1,2,k))        &
            +tr12*(cc(m,i-1,5,k)+cc(m,ic-1,4,k)))+(ti11*(cc(m,i,3,k)    &
            +cc(m,ic,2,k))+ti12*(cc(m,i,5,k)+cc(m,ic,4,k))))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE vradb5


SUBROUTINE vradbg (mp,ido,ip,l1,idl1,cc,c1,c2,ch,ch2,                   & 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
! &
         mdimc,wa)
  DIMENSION    ch(mdimc,ido,l1,ip)    ,cc(mdimc,ido,ip,l1) ,            &
             c1(mdimc,ido,l1,ip)     ,c2(mdimc,idl1,ip),                &
                  ch2(mdimc,idl1,ip)       ,wa(ido)
  tpi=2.*pimach(1.0)
  arg = tpi/FLOAT(ip)
  dcp = COS(arg)
  dsp = SIN(arg)
  idp2 = ido+2
  nbd = (ido-1)/2
  ipp2 = ip+2
  ipph = (ip+1)/2
  IF (ido < l1) GO TO 103
  DO k=1,l1
    DO i=1,ido
      DO m=1,mp
        ch(m,i,k,1) = cc(m,i,1,k)
      END DO
    END DO
  END DO
  GO TO 106
  103 DO i=1,ido
    DO k=1,l1
      DO m=1,mp
        ch(m,i,k,1) = cc(m,i,1,k)
      END DO
    END DO
  END DO
  106 DO j=2,ipph
    jc = ipp2-j
    j2 = j+j
    DO k=1,l1
      DO m=1,mp
        ch(m,1,k,j) = cc(m,ido,j2-2,k)+cc(m,ido,j2-2,k)
        ch(m,1,k,jc) = cc(m,1,j2-1,k)+cc(m,1,j2-1,k)
      END DO
    END DO
  END DO
  IF (ido == 1) GO TO 116
  IF (nbd < l1) GO TO 112
  DO j=2,ipph
    jc = ipp2-j
    DO k=1,l1
      DO i=3,ido,2
        ic = idp2-i
        DO m=1,mp
          ch(m,i-1,k,j) = cc(m,i-1,2*j-1,k)+cc(m,ic-1,2*j-2,k)
          ch(m,i-1,k,jc) = cc(m,i-1,2*j-1,k)-cc(m,ic-1,2*j-2,k)
          ch(m,i,k,j) = cc(m,i,2*j-1,k)-cc(m,ic,2*j-2,k)
          ch(m,i,k,jc) = cc(m,i,2*j-1,k)+cc(m,ic,2*j-2,k)
        END DO
      END DO
    END DO
  END DO
  GO TO 116
  112 DO j=2,ipph
    jc = ipp2-j
    DO i=3,ido,2
      ic = idp2-i
      DO k=1,l1
        DO m=1,mp
          ch(m,i-1,k,j) = cc(m,i-1,2*j-1,k)+cc(m,ic-1,2*j-2,k)
          ch(m,i-1,k,jc) = cc(m,i-1,2*j-1,k)-cc(m,ic-1,2*j-2,k)
          ch(m,i,k,j) = cc(m,i,2*j-1,k)-cc(m,ic,2*j-2,k)
          ch(m,i,k,jc) = cc(m,i,2*j-1,k)+cc(m,ic,2*j-2,k)
        END DO
      END DO
    END DO
  END DO
  116 ar1 = 1.
  ai1 = 0.
  DO l=2,ipph
    lc = ipp2-l
    ar1h = dcp*ar1-dsp*ai1
    ai1 = dcp*ai1+dsp*ar1
    ar1 = ar1h
    DO ik=1,idl1
      DO m=1,mp
        c2(m,ik,l) = ch2(m,ik,1)+ar1*ch2(m,ik,2)
        c2(m,ik,lc) = ai1*ch2(m,ik,ip)
      END DO
    END DO
    dc2 = ar1
    ds2 = ai1
    ar2 = ar1
    ai2 = ai1
    DO j=3,ipph
      jc = ipp2-j
      ar2h = dc2*ar2-ds2*ai2
      ai2 = dc2*ai2+ds2*ar2
      ar2 = ar2h
      DO ik=1,idl1
        DO m=1,mp
          c2(m,ik,l) = c2(m,ik,l)+ar2*ch2(m,ik,j)
          c2(m,ik,lc) = c2(m,ik,lc)+ai2*ch2(m,ik,jc)
        END DO
      END DO
    END DO
  END DO
  DO j=2,ipph
    DO ik=1,idl1
      DO m=1,mp
        ch2(m,ik,1) = ch2(m,ik,1)+ch2(m,ik,j)
      END DO
    END DO
  END DO
  DO j=2,ipph
    jc = ipp2-j
    DO k=1,l1
      DO m=1,mp
        ch(m,1,k,j) = c1(m,1,k,j)-c1(m,1,k,jc)
        ch(m,1,k,jc) = c1(m,1,k,j)+c1(m,1,k,jc)
      END DO
    END DO
  END DO
  IF (ido == 1) GO TO 132
  IF (nbd < l1) GO TO 128
  DO j=2,ipph
    jc = ipp2-j
    DO k=1,l1
      DO i=3,ido,2
        DO m=1,mp
          ch(m,i-1,k,j) = c1(m,i-1,k,j)-c1(m,i,k,jc)
          ch(m,i-1,k,jc) = c1(m,i-1,k,j)+c1(m,i,k,jc)
          ch(m,i,k,j) = c1(m,i,k,j)+c1(m,i-1,k,jc)
          ch(m,i,k,jc) = c1(m,i,k,j)-c1(m,i-1,k,jc)
        END DO
      END DO
    END DO
  END DO
  GO TO 132
  128 DO j=2,ipph
    jc = ipp2-j
    DO i=3,ido,2
      DO k=1,l1
        DO m=1,mp
          ch(m,i-1,k,j) = c1(m,i-1,k,j)-c1(m,i,k,jc)
          ch(m,i-1,k,jc) = c1(m,i-1,k,j)+c1(m,i,k,jc)
          ch(m,i,k,j) = c1(m,i,k,j)+c1(m,i-1,k,jc)
          ch(m,i,k,jc) = c1(m,i,k,j)-c1(m,i-1,k,jc)
        END DO
      END DO
    END DO
  END DO
  132 CONTINUE
  IF (ido == 1) RETURN
  DO ik=1,idl1
    DO m=1,mp
      c2(m,ik,1) = ch2(m,ik,1)
    END DO
  END DO
  DO j=2,ip
    DO k=1,l1
      DO m=1,mp
        c1(m,1,k,j) = ch(m,1,k,j)
      END DO
    END DO
  END DO
  IF (nbd > l1) GO TO 139
  is = -ido
  DO j=2,ip
    is = is+ido
    idij = is
    DO i=3,ido,2
      idij = idij+2
      DO k=1,l1
        DO m=1,mp
          c1(m,i-1,k,j) = wa(idij-1)*ch(m,i-1,k,j)-wa(idij)*            &
              ch(m,i,k,j)
          c1(m,i,k,j) = wa(idij-1)*ch(m,i,k,j)+wa(idij)*                &
              ch(m,i-1,k,j)
        END DO
      END DO
    END DO
  END DO
  GO TO 143
  139 is = -ido
  DO j=2,ip
    is = is+ido
    DO k=1,l1
      idij = is
      DO i=3,ido,2
        idij = idij+2
        DO m=1,mp
          c1(m,i-1,k,j) = wa(idij-1)*ch(m,i-1,k,j)-wa(idij)*            &
              ch(m,i,k,j)
          c1(m,i,k,j) = wa(idij-1)*ch(m,i,k,j)+wa(idij)*                &
              ch(m,i-1,k,j)
        END DO
      END DO
    END DO
  END DO
  143 RETURN
END SUBROUTINE vradbg


SUBROUTINE vradf2 (mp,ido,l1,cc,ch,mdimc,wa1) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION   ch(mdimc,ido,2,l1)  ,cc(mdimc,ido,l1,2)     ,             &
                  wa1(ido)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,1,k) = cc(m,1,k,1)+cc(m,1,k,2)
      ch(m,ido,2,k) = cc(m,1,k,1)-cc(m,1,k,2)
    END DO
  END DO
  IF (ido-2 < 0) THEN
    GO TO   107
  ELSE IF (ido-2 == 0) THEN
    GO TO   105
  END IF
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i,1,k) = cc(m,i,k,1)+(wa1(i-2)*cc(m,i,k,2)-                &
            wa1(i-1)*cc(m,i-1,k,2))
        ch(m,ic,2,k) = (wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*                  &
            cc(m,i-1,k,2))-cc(m,i,k,1)
        ch(m,i-1,1,k) = cc(m,i-1,k,1)+(wa1(i-2)*cc(m,i-1,k,2)+          &
            wa1(i-1)*cc(m,i,k,2))
        ch(m,ic-1,2,k) = cc(m,i-1,k,1)-(wa1(i-2)*cc(m,i-1,k,2)+         &
            wa1(i-1)*cc(m,i,k,2))
      END DO
    END DO
  END DO
  IF (MOD(ido,2) == 1) RETURN
  105 DO k=1,l1
    DO m=1,mp
      ch(m,1,2,k) = -cc(m,ido,k,2)
      ch(m,ido,1,k) = cc(m,ido,k,1)
    END DO
  END DO
  107 RETURN
END SUBROUTINE vradf2


SUBROUTINE vradf3 (mp,ido,l1,cc,ch,mdimc,wa1,wa2) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION   ch(mdimc,ido,3,l1)  ,cc(mdimc,ido,l1,3)     ,             &
                  wa1(ido)     ,wa2(ido)
  arg=2.*pimach(1.0)/3.
  taur=COS(arg)
  taui=SIN(arg)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,1,k) = cc(m,1,k,1)+(cc(m,1,k,2)+cc(m,1,k,3))
      ch(m,1,3,k) = taui*(cc(m,1,k,3)-cc(m,1,k,2))
      ch(m,ido,2,k) = cc(m,1,k,1)+taur*                                 &
          (cc(m,1,k,2)+cc(m,1,k,3))
    END DO
  END DO
  IF (ido == 1) RETURN
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,1,k) = cc(m,i-1,k,1)+((wa1(i-2)*cc(m,i-1,k,2)+         &
            wa1(i-1)*cc(m,i,k,2))+(wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*     &
            cc(m,i,k,3)))
        ch(m,i,1,k) = cc(m,i,k,1)+((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*      &
            cc(m,i-1,k,2))+(wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*              &
            cc(m,i-1,k,3)))
        ch(m,i-1,3,k) = (cc(m,i-1,k,1)+taur*((wa1(i-2)*                 &
            cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,k,2))+(wa2(i-2)*              &
            cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,k,3))))+(taui*((wa1(i-2)*     &
            cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,2))-(wa2(i-2)*              &
            cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k,3))))
        ch(m,ic-1,2,k) = (cc(m,i-1,k,1)+taur*((wa1(i-2)*                &
            cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,k,2))+(wa2(i-2)*              &
            cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,k,3))))-(taui*((wa1(i-2)*     &
            cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,2))-(wa2(i-2)*              &
            cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k,3))))
        ch(m,i,3,k) = (cc(m,i,k,1)+taur*((wa1(i-2)*cc(m,i,k,2)-         &
            wa1(i-1)*cc(m,i-1,k,2))+(wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*     &
            cc(m,i-1,k,3))))+(taui*((wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*   &
            cc(m,i,k,3))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2))))
        ch(m,ic,2,k) = (taui*((wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*         &
            cc(m,i,k,3))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2))))-(cc(m,i,k,1)+taur*((wa1(i-2)*cc(m,i,k,2)-    &
            wa1(i-1)*cc(m,i-1,k,2))+(wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*     &
            cc(m,i-1,k,3))))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE vradf3


SUBROUTINE vradf4 (mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION    cc(mdimc,ido,l1,4)   ,ch(mdimc,ido,4,l1)     ,           &
                  wa1(ido)     ,wa2(ido)     ,wa3(ido)
  hsqt2=SQRT(2.)/2.
  DO k=1,l1
    DO m=1,mp
      ch(m,1,1,k) = (cc(m,1,k,2)+cc(m,1,k,4))                           &
          +(cc(m,1,k,1)+cc(m,1,k,3))
      ch(m,ido,4,k) = (cc(m,1,k,1)+cc(m,1,k,3))                         &
          -(cc(m,1,k,2)+cc(m,1,k,4))
      ch(m,ido,2,k) = cc(m,1,k,1)-cc(m,1,k,3)
      ch(m,1,3,k) = cc(m,1,k,4)-cc(m,1,k,2)
    END DO
  END DO
  IF (ido-2 < 0) THEN
    GO TO   107
  ELSE IF (ido-2 == 0) THEN
    GO TO   105
  END IF
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,1,k) = ((wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2))+(wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*              &
            cc(m,i,k,4)))+(cc(m,i-1,k,1)+(wa2(i-2)*cc(m,i-1,k,3)+       &
            wa2(i-1)*cc(m,i,k,3)))
        ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+(wa2(i-2)*cc(m,i-1,k,3)+        &
            wa2(i-1)*cc(m,i,k,3)))-((wa1(i-2)*cc(m,i-1,k,2)+            &
            wa1(i-1)*cc(m,i,k,2))+(wa3(i-2)*cc(m,i-1,k,4)+              &
            wa3(i-1)*cc(m,i,k,4)))
        ch(m,i,1,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*                  &
            cc(m,i-1,k,2))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4)))+(cc(m,i,k,1)+(wa2(i-2)*cc(m,i,k,3)-         &
            wa2(i-1)*cc(m,i-1,k,3)))
        ch(m,ic,4,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*                 &
            cc(m,i-1,k,2))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4)))-(cc(m,i,k,1)+(wa2(i-2)*cc(m,i,k,3)-         &
            wa2(i-1)*cc(m,i-1,k,3)))
        ch(m,i-1,3,k) = ((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*                &
            cc(m,i-1,k,2))-(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4)))+(cc(m,i-1,k,1)-(wa2(i-2)*cc(m,i-1,k,3)+     &
            wa2(i-1)*cc(m,i,k,3)))
        ch(m,ic-1,2,k) = (cc(m,i-1,k,1)-(wa2(i-2)*cc(m,i-1,k,3)+        &
            wa2(i-1)*cc(m,i,k,3)))-((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*     &
            cc(m,i-1,k,2))-(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4)))
        ch(m,i,3,k) = ((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*                &
            cc(m,i,k,4))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2)))+(cc(m,i,k,1)-(wa2(i-2)*cc(m,i,k,3)-           &
            wa2(i-1)*cc(m,i-1,k,3)))
        ch(m,ic,2,k) = ((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*               &
            cc(m,i,k,4))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2)))-(cc(m,i,k,1)-(wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*  &
            cc(m,i-1,k,3)))
      END DO
    END DO
  END DO
  IF (MOD(ido,2) == 1) RETURN
  105 CONTINUE
  DO k=1,l1
    DO m=1,mp
      ch(m,ido,1,k) = (hsqt2*(cc(m,ido,k,2)-cc(m,ido,k,4)))+            &
          cc(m,ido,k,1)
      ch(m,ido,3,k) = cc(m,ido,k,1)-(hsqt2*(cc(m,ido,k,2)-              &
          cc(m,ido,k,4)))
      ch(m,1,2,k) = (-hsqt2*(cc(m,ido,k,2)+cc(m,ido,k,4)))-             &
          cc(m,ido,k,3)
      ch(m,1,4,k) = (-hsqt2*(cc(m,ido,k,2)+cc(m,ido,k,4)))+             &
          cc(m,ido,k,3)
    END DO
  END DO
  107 RETURN
END SUBROUTINE vradf4


SUBROUTINE vradf5 (mp,ido,l1,cc,ch,mdimc,wa1,wa2,wa3,wa4) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION  cc(mdimc,ido,l1,5)    ,ch(mdimc,ido,5,l1)     ,            &
             wa1(ido)     ,wa2(ido)     ,wa3(ido)     ,wa4(ido)
  arg=2.*pimach(1.0)/5.
  tr11=COS(arg)
  ti11=SIN(arg)
  tr12=COS(2.*arg)
  ti12=SIN(2.*arg)
  DO k=1,l1
    DO m=1,mp
      ch(m,1,1,k) = cc(m,1,k,1)+(cc(m,1,k,5)+cc(m,1,k,2))+              &
          (cc(m,1,k,4)+cc(m,1,k,3))
      ch(m,ido,2,k) = cc(m,1,k,1)+tr11*(cc(m,1,k,5)+cc(m,1,k,2))+       &
          tr12*(cc(m,1,k,4)+cc(m,1,k,3))
      ch(m,1,3,k) = ti11*(cc(m,1,k,5)-cc(m,1,k,2))+ti12*                &
          (cc(m,1,k,4)-cc(m,1,k,3))
      ch(m,ido,4,k) = cc(m,1,k,1)+tr12*(cc(m,1,k,5)+cc(m,1,k,2))+       &
          tr11*(cc(m,1,k,4)+cc(m,1,k,3))
      ch(m,1,5,k) = ti12*(cc(m,1,k,5)-cc(m,1,k,2))-ti11*                &
          (cc(m,1,k,4)-cc(m,1,k,3))
    END DO
  END DO
  IF (ido == 1) RETURN
  idp2 = ido+2
  DO k=1,l1
    DO i=3,ido,2
      ic = idp2-i
      DO m=1,mp
        ch(m,i-1,1,k) = cc(m,i-1,k,1)+((wa1(i-2)*cc(m,i-1,k,2)+         &
            wa1(i-1)*cc(m,i,k,2))+(wa4(i-2)*cc(m,i-1,k,5)+wa4(i-1)*     &
            cc(m,i,k,5)))+((wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*            &
            cc(m,i,k,3))+(wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k,4)))
        ch(m,i,1,k) = cc(m,i,k,1)+((wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*      &
            cc(m,i-1,k,2))+(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*              &
            cc(m,i-1,k,5)))+((wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*            &
            cc(m,i-1,k,3))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4)))
        ch(m,i-1,3,k) = cc(m,i-1,k,1)+tr11*                             &
            ( wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,k,2)               &
            +wa4(i-2)*cc(m,i-1,k,5)+wa4(i-1)*cc(m,i,k,5))+tr12*         &
            ( wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,k,3)               &
            +wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k,4))+ti11*         &
            ( wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,2)               &
            -(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m,i-1,k,5)))+ti12*       &
            ( wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k,3)               &
            -(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1,k,4)))
        ch(m,ic-1,2,k) = cc(m,i-1,k,1)+tr11*                            &
            ( wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,k,2)               &
            +wa4(i-2)*cc(m,i-1,k,5)+wa4(i-1)*cc(m,i,k,5))+tr12*         &
            ( wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,k,3)               &
            +wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k,4))-(ti11*        &
            ( wa1(i-2)*cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,2)               &
            -(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*cc(m,i-1,k,5)))+ti12*       &
            ( wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*cc(m,i-1,k,3)               &
            -(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*cc(m,i-1,k,4))))
        ch(m,i,3,k) = (cc(m,i,k,1)+tr11*((wa1(i-2)*cc(m,i,k,2)-         &
            wa1(i-1)*cc(m,i-1,k,2))+(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*     &
            cc(m,i-1,k,5)))+tr12*((wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*       &
            cc(m,i-1,k,3))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4))))+(ti11*((wa4(i-2)*cc(m,i-1,k,5)+            &
            wa4(i-1)*cc(m,i,k,5))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*     &
            cc(m,i,k,2)))+ti12*((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*       &
            cc(m,i,k,4))-(wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*              &
            cc(m,i,k,3))))
        ch(m,ic,2,k) = (ti11*((wa4(i-2)*cc(m,i-1,k,5)+wa4(i-1)*         &
            cc(m,i,k,5))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2)))+ti12*((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*       &
            cc(m,i,k,4))-(wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*              &
            cc(m,i,k,3))))-(cc(m,i,k,1)+tr11*((wa1(i-2)*cc(m,i,k,2)-    &
            wa1(i-1)*cc(m,i-1,k,2))+(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*     &
            cc(m,i-1,k,5)))+tr12*((wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*       &
            cc(m,i-1,k,3))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4))))
        ch(m,i-1,5,k) = (cc(m,i-1,k,1)+tr12*((wa1(i-2)*                 &
            cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,k,2))+(wa4(i-2)*              &
            cc(m,i-1,k,5)+wa4(i-1)*cc(m,i,k,5)))+tr11*((wa2(i-2)*       &
            cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,k,3))+(wa3(i-2)*              &
            cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k,4))))+(ti12*((wa1(i-2)*     &
            cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,2))-(wa4(i-2)*cc(m,i,k,5)-  &
            wa4(i-1)*cc(m,i-1,k,5)))-ti11*((wa2(i-2)*cc(m,i,k,3)-       &
            wa2(i-1)*cc(m,i-1,k,3))-(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*     &
            cc(m,i-1,k,4))))
        ch(m,ic-1,4,k) = (cc(m,i-1,k,1)+tr12*((wa1(i-2)*                &
            cc(m,i-1,k,2)+wa1(i-1)*cc(m,i,k,2))+(wa4(i-2)*              &
            cc(m,i-1,k,5)+wa4(i-1)*cc(m,i,k,5)))+tr11*((wa2(i-2)*       &
            cc(m,i-1,k,3)+wa2(i-1)*cc(m,i,k,3))+(wa3(i-2)*              &
            cc(m,i-1,k,4)+wa3(i-1)*cc(m,i,k,4))))-(ti12*((wa1(i-2)*     &
            cc(m,i,k,2)-wa1(i-1)*cc(m,i-1,k,2))-(wa4(i-2)*cc(m,i,k,5)-  &
            wa4(i-1)*cc(m,i-1,k,5)))-ti11*((wa2(i-2)*cc(m,i,k,3)-       &
            wa2(i-1)*cc(m,i-1,k,3))-(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*     &
            cc(m,i-1,k,4))))
        ch(m,i,5,k) = (cc(m,i,k,1)+tr12*((wa1(i-2)*cc(m,i,k,2)-         &
            wa1(i-1)*cc(m,i-1,k,2))+(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*     &
            cc(m,i-1,k,5)))+tr11*((wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*       &
            cc(m,i-1,k,3))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4))))+(ti12*((wa4(i-2)*cc(m,i-1,k,5)+            &
            wa4(i-1)*cc(m,i,k,5))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*     &
            cc(m,i,k,2)))-ti11*((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*       &
            cc(m,i,k,4))-(wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*              &
            cc(m,i,k,3))))
        ch(m,ic,4,k) = (ti12*((wa4(i-2)*cc(m,i-1,k,5)+wa4(i-1)*         &
            cc(m,i,k,5))-(wa1(i-2)*cc(m,i-1,k,2)+wa1(i-1)*              &
            cc(m,i,k,2)))-ti11*((wa3(i-2)*cc(m,i-1,k,4)+wa3(i-1)*       &
            cc(m,i,k,4))-(wa2(i-2)*cc(m,i-1,k,3)+wa2(i-1)*              &
            cc(m,i,k,3))))-(cc(m,i,k,1)+tr12*((wa1(i-2)*cc(m,i,k,2)-    &
            wa1(i-1)*cc(m,i-1,k,2))+(wa4(i-2)*cc(m,i,k,5)-wa4(i-1)*     &
            cc(m,i-1,k,5)))+tr11*((wa2(i-2)*cc(m,i,k,3)-wa2(i-1)*       &
            cc(m,i-1,k,3))+(wa3(i-2)*cc(m,i,k,4)-wa3(i-1)*              &
            cc(m,i-1,k,4))))
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE vradf5


SUBROUTINE vradfg (mp,ido,ip,l1,idl1,cc,c1,c2,ch,ch2,mdimc,wa) 2
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION     ch(mdimc,ido,l1,ip)   ,cc(mdimc,ido,ip,l1)  ,           &
              c1(mdimc,ido,l1,ip)    ,c2(mdimc,idl1,ip),                &
                  ch2(mdimc,idl1,ip)           ,wa(ido)
  tpi=2.*pimach(1.0)
  arg = tpi/FLOAT(ip)
  dcp = COS(arg)
  dsp = SIN(arg)
  ipph = (ip+1)/2
  ipp2 = ip+2
  idp2 = ido+2
  nbd = (ido-1)/2
  IF (ido == 1) GO TO 119
  DO ik=1,idl1
    DO m=1,mp
      ch2(m,ik,1) = c2(m,ik,1)
    END DO
  END DO
  DO j=2,ip
    DO k=1,l1
      DO m=1,mp
        ch(m,1,k,j) = c1(m,1,k,j)
      END DO
    END DO
  END DO
  IF (nbd > l1) GO TO 107
  is = -ido
  DO j=2,ip
    is = is+ido
    idij = is
    DO i=3,ido,2
      idij = idij+2
      DO k=1,l1
        DO m=1,mp
          ch(m,i-1,k,j) = wa(idij-1)*c1(m,i-1,k,j)+wa(idij)             &
              *c1(m,i,k,j)
          ch(m,i,k,j) = wa(idij-1)*c1(m,i,k,j)-wa(idij)                 &
              *c1(m,i-1,k,j)
        END DO
      END DO
    END DO
  END DO
  GO TO 111
  107 is = -ido
  DO j=2,ip
    is = is+ido
    DO k=1,l1
      idij = is
      DO i=3,ido,2
        idij = idij+2
        DO m=1,mp
          ch(m,i-1,k,j) = wa(idij-1)*c1(m,i-1,k,j)+wa(idij)             &
              *c1(m,i,k,j)
          ch(m,i,k,j) = wa(idij-1)*c1(m,i,k,j)-wa(idij)                 &
              *c1(m,i-1,k,j)
        END DO
      END DO
    END DO
  END DO
  111 IF (nbd < l1) GO TO 115
  DO j=2,ipph
    jc = ipp2-j
    DO k=1,l1
      DO i=3,ido,2
        DO m=1,mp
          c1(m,i-1,k,j) = ch(m,i-1,k,j)+ch(m,i-1,k,jc)
          c1(m,i-1,k,jc) = ch(m,i,k,j)-ch(m,i,k,jc)
          c1(m,i,k,j) = ch(m,i,k,j)+ch(m,i,k,jc)
          c1(m,i,k,jc) = ch(m,i-1,k,jc)-ch(m,i-1,k,j)
        END DO
      END DO
    END DO
  END DO
  GO TO 121
  115 DO j=2,ipph
    jc = ipp2-j
    DO i=3,ido,2
      DO k=1,l1
        DO m=1,mp
          c1(m,i-1,k,j) = ch(m,i-1,k,j)+ch(m,i-1,k,jc)
          c1(m,i-1,k,jc) = ch(m,i,k,j)-ch(m,i,k,jc)
          c1(m,i,k,j) = ch(m,i,k,j)+ch(m,i,k,jc)
          c1(m,i,k,jc) = ch(m,i-1,k,jc)-ch(m,i-1,k,j)
        END DO
      END DO
    END DO
  END DO
  GO TO 121
  119 DO ik=1,idl1
    DO m=1,mp
      c2(m,ik,1) = ch2(m,ik,1)
    END DO
  END DO
  121 DO j=2,ipph
    jc = ipp2-j
    DO k=1,l1
      DO m=1,mp
        c1(m,1,k,j) = ch(m,1,k,j)+ch(m,1,k,jc)
        c1(m,1,k,jc) = ch(m,1,k,jc)-ch(m,1,k,j)
      END DO
    END DO
  END DO
!
  ar1 = 1.
  ai1 = 0.
  DO l=2,ipph
    lc = ipp2-l
    ar1h = dcp*ar1-dsp*ai1
    ai1 = dcp*ai1+dsp*ar1
    ar1 = ar1h
    DO ik=1,idl1
      DO m=1,mp
        ch2(m,ik,l) = c2(m,ik,1)+ar1*c2(m,ik,2)
        ch2(m,ik,lc) = ai1*c2(m,ik,ip)
      END DO
    END DO
    dc2 = ar1
    ds2 = ai1
    ar2 = ar1
    ai2 = ai1
    DO j=3,ipph
      jc = ipp2-j
      ar2h = dc2*ar2-ds2*ai2
      ai2 = dc2*ai2+ds2*ar2
      ar2 = ar2h
      DO ik=1,idl1
        DO m=1,mp
          ch2(m,ik,l) = ch2(m,ik,l)+ar2*c2(m,ik,j)
          ch2(m,ik,lc) = ch2(m,ik,lc)+ai2*c2(m,ik,jc)
        END DO
      END DO
    END DO
  END DO
  DO j=2,ipph
    DO ik=1,idl1
      DO m=1,mp
        ch2(m,ik,1) = ch2(m,ik,1)+c2(m,ik,j)
      END DO
    END DO
  END DO
!
  IF (ido < l1) GO TO 132
  DO k=1,l1
    DO i=1,ido
      DO m=1,mp
        cc(m,i,1,k) = ch(m,i,k,1)
      END DO
    END DO
  END DO
  GO TO 135
  132 DO i=1,ido
    DO k=1,l1
      DO m=1,mp
        cc(m,i,1,k) = ch(m,i,k,1)
      END DO
    END DO
  END DO
  135 DO j=2,ipph
    jc = ipp2-j
    j2 = j+j
    DO k=1,l1
      DO m=1,mp
        cc(m,ido,j2-2,k) = ch(m,1,k,j)
        cc(m,1,j2-1,k) = ch(m,1,k,jc)
      END DO
    END DO
  END DO
  IF (ido == 1) RETURN
  IF (nbd < l1) GO TO 141
  DO j=2,ipph
    jc = ipp2-j
    j2 = j+j
    DO k=1,l1
      DO i=3,ido,2
        ic = idp2-i
        DO m=1,mp
          cc(m,i-1,j2-1,k) = ch(m,i-1,k,j)+ch(m,i-1,k,jc)
          cc(m,ic-1,j2-2,k) = ch(m,i-1,k,j)-ch(m,i-1,k,jc)
          cc(m,i,j2-1,k) = ch(m,i,k,j)+ch(m,i,k,jc)
          cc(m,ic,j2-2,k) = ch(m,i,k,jc)-ch(m,i,k,j)
        END DO
      END DO
    END DO
  END DO
  RETURN
  141 DO j=2,ipph
    jc = ipp2-j
    j2 = j+j
    DO i=3,ido,2
      ic = idp2-i
      DO k=1,l1
        DO m=1,mp
          cc(m,i-1,j2-1,k) = ch(m,i-1,k,j)+ch(m,i-1,k,jc)
          cc(m,ic-1,j2-2,k) = ch(m,i-1,k,j)-ch(m,i-1,k,jc)
          cc(m,i,j2-1,k) = ch(m,i,k,j)+ch(m,i,k,jc)
          cc(m,ic,j2-2,k) = ch(m,i,k,jc)-ch(m,i,k,j)
        END DO
      END DO
    END DO
  END DO
  RETURN
END SUBROUTINE vradfg


SUBROUTINE vrfftb(m,n,r,rt,mdimr,wsave),1
!***BEGIN PROLOGUE  VRFFTB
!***DATE WRITTEN   850801   (YYMMDD)
!***REVISION DATE  900509   (YYMMDD)
!***CATEGORY NO.  J1A1
!***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
!          FOURIER SYNTHESIS, BACKWARD TRANSFORM, MULTIPLE SEQUENCES
!***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
!***PURPOSE  Backward real periodic transform, M sequences.
!***DESCRIPTION
!
!  Subroutine VRFFTB computes the synthesis (backward transform) of a
!  number of real periodic sequences from their Fourier coefficients.
!  Specifically, for each set of independent Fourier coefficients
!  F(K), the corresponding real periodic sequence is computed.
!
!  The array WSAVE which is used by subroutine VRFFTB must be
!  initialized by calling subroutine VRFFTI(N,WSAVE).
!
!
!  Input Parameters
!
!  M       the number of sets of coefficients.
!
!  N       the length of the sequences of coefficients to be
!       transformed.  The method is most efficient when N is a
!       product of small primes, however n may be any positive
!       integer.
!
!  R       areal two-dimensional array of size MDIMX x N containing the
!       coefficients to be transformed.  Each set of coefficients
!       F(K), K\0,1,..,N-1, is stored as a ROW of R.  Specifically,
!       the I-th set of independent Fourier coefficients is stored
!
!             R(I,1) = REAL( F(I,0) ),
!
!             R(I,2*K) = REAL( F(I,K) )
!
!             R(I,2*K+1) = IMAG( F(I,K) )
!
!                for K = 1, 2, . . . , M-1,
!
!             and, when N is even,
!
!             R(I,N) = REAL( F(I,N/2) ).
!
!  RT      a real two-dimensional work array of size MDIMX x N.
!
!  MDIMR   the row (or first) dimension of the arrays R and RT exactly
!       as they appear in the calling program.  This parameter is
!       used to specify the variable dimension of these arrays.
!
!  WSAVE   a real one-dimensional work array which must be dimensioned
!       at least N+15.  The WSAVE array must be initialized by
!       calling subroutine VRFFTI.  A different WSAVE array must be
!       used for each different value of N.  This initialization does
!       not have to be repeated so long as N remains unchanged.  The
!       same WSAVE array may be used by VRFFTB and VRFFTB.
!
!  Output Parameters
!
!  R       contains M real periodic sequences corresponding to the given
!       coefficients.  Specifically, the I-th row of R contains the
!       real periodic sequence corresponding to the I-th set of
!       independent Fourier coefficients F(I,K) stored as
!
!            R(I,J) = X(I,J-1) ,   J = 1, 2, . . . , N, where
!
!            X(I,J) = SQRT(1/N)* F(I,0) + (-1)**J*F(I,N/2)
!                     + 2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N)
!                     - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ]  ,
!
!              when N is even, and
!
!            X(I,J) = SQRT(1/N)* F(I,0) +
!                     2*SUM(K=1,M)[ REAL(F(I,2K))*COS(2K*J*PI/N)
!                     - IMAG(F(I,2K+1))*SIN(2K*J*PI/N) ]  ,
!
!              when N is odd.
!
!  WSAVE   contains results which must not be destroyed between calls
!       to VRFFTF or VRFFTB.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VRFFTF followed immediately by a call of
!        of VRFFTB will return the original sequences R.  Thus,
!        VRFFTB is the correctly normalized inverse of VRFFTF.
!
!  -----------------------------------------------------------------
!
!  VRFFTB is a straightforward extension of the subprogram RFFTB to
!  handle M simultaneous sequences.  RFFTB was originally developed
!  by P. N. Swarztrauber of NCAR.
!
!
!           * * * * * * * * * * * * * * * * * * * * *
!           *                                       *
!           *         PROGRAM SPECIFICATIONS        *
!           *                                       *
!           * * * * * * * * * * * * * * * * * * * * *
!
!
!  DIMENSION OF    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
!  ARGUMENTS
!
!  LATEST          AUGUST 1, 1985
!  REVISION
!
!  SUBPROGRAMS     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
!  REQUIRED        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
!                  VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
!
!  SPECIAL         NONE
!  CONDITIONS
!
!  COMMON          NONE
!  BLOCKS
!
!  I/O             NONE
!
!  PRECISION       SINGLE
!
!  SPECIALIST      ROLAND SWEET
!
!  LANGUAGE        FORTRAN
!
!  HISTORY         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
!                  NATIONAL BUREAU OF STANDARDS (BOULDER).
!
!  ALGORITHM       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
!                  OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
!
!  PORTABILITY     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
!                  THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
!                  THE FUNCTION PIMACH.
!
!  REQUIRED        COS,SIN
!  RESIDENT
!  ROUTINES
!
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!            Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!            pp. 51-83.
!***ROUTINES CALLED  VRFTB1
!***END PROLOGUE  VRFFTB
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION     r(mdimr,n),rt(mdimr,n),wsave(n+15)
  IF (n == 1) RETURN
  CALL vrftb1 (m,n,r,rt,mdimr,wsave(1),wsave(n+1))
  RETURN
END SUBROUTINE vrfftb


SUBROUTINE vrfftf (m,n,r,rt,mdimr,wsave) 1,1
!***BEGIN PROLOGUE  VRFFTF
!***DATE WRITTEN   850801   (YYMMDD)
!***REVISION DATE  900509   (YYMMDD)
!***CATEGORY NO.  J1A1
!***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
!          FOURIER ANALYSIS, FORWARD TRANSFORM, MULTIPLE SEQUENCES
!***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
!***PURPOSE  Forward real periodic transform, M sequences.
!***DESCRIPTION
!
!  Subroutine VRFFTF computes the Fourier coefficients (forward
!  transform) of a number of real periodic sequences.  Specifically,
!  for each sequence the subroutine claculates the independent
!  Fourier coefficients described below at output parameter R.
!
!  The array WSAVE which is used by subroutine VRFFTF must be
!  initialized by calling subroutine VRFFTI(N,WSAVE).
!
!
!  Input Parameters
!
!  M       the number of sequences to be transformed.
!
!  N       the length of the sequences to be transformed.  The method
!       is most efficient when N is a product of small primes,
!       however n may be any positive integer.
!
!  R       areal two-dimensional array of size MDIMX x N containing the
!       the sequences to be transformed.  The sequences are stored
!       in the ROWS of R.  Thus, the I-th sequence to be transformed,
!       X(I,J), J=0,1,...,N-1, is stored as
!
!            R(I,J) = X(I,J-1) , J=1, 2, . . . , N.
!
!  RT      a real two-dimensional work array of size MDIMX x N.
!
!  MDIMR   the row (or first) dimension of the arrays R and RT exactly
!       as they appear in the calling program.  This parameter is
!       used to specify the variable dimension of these arrays.
!
!  WSAVE   a real one-dimensional work array which must be dimensioned
!       at least N+15.  The WSAVE array must be initialized by
!       calling subroutine VRFFTI.  A different WSAVE array must be
!       used for each different value of N.  This initialization does
!       not have to be repeated so long as N remains unchanged.  The
!       same WSAVE array may be used by VRFFTF and VRFFTB.
!
!  Output Parameters
!
!  R       contains the Fourier coefficients F(K) for each of the M
!       input sequences.  Specifically, row I of R, R(I,J),
!       J=1,2,..,N, contains the independent Fourier coefficients
!       F(I,K), for the I-th input sequence stored as
!
!          R(I,1) = REAL( F(I,0) ),
!                 = SQRT(1/N)*SUM(J=0,N-1)[ X(I,J) ],
!
!          R(I,2*K) = REAL( F(I,K) )
!                   = SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*COS(2J*K*PI/N)]
!
!          R(I,2*K+1) = IMAG( F(I,K) )
!                     =-SQRT(1/N)*SUM(J=0,N-1)[X(I,J)*SIN(2J*K*PI/N)]
!
!                for K = 1, 2, . . . , M-1,
!
!           and, when N is even,
!
!           R(I,N) = REAL( F(I,N/2) ).
!                  = SQRT(1/N)*SUM(J=0,N-1)[ (-1)**J*X(I,J) ].
!
!  WSAVE   contains results which must not be destroyed between calls
!       to VRFFTF or VRFFTB.
!
!  -----------------------------------------------------------------
!
!  NOTE  -  A call of VRFFTF followed immediately by a call of
!        of VRFFTB will return the original sequences R.  Thus,
!        VRFFTB is the correctly normalized inverse of VRFFTF.
!
!  -----------------------------------------------------------------
!
!  VRFFTF is a straightforward extension of the subprogram RFFTF to
!  handle M simultaneous sequences.  RFFTF was originally developed
!  by P. N. Swarztrauber of NCAR.
!
!
!           * * * * * * * * * * * * * * * * * * * * *
!           *                                       *
!           *         PROGRAM SPECIFICATIONS        *
!           *                                       *
!           * * * * * * * * * * * * * * * * * * * * *
!
!
!  DIMENSION OF    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
!  ARGUMENTS
!
!  LATEST          AUGUST 1, 1985
!  REVISION
!
!  SUBPROGRAMS     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
!  REQUIRED        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
!                  VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
!
!  SPECIAL         NONE
!  CONDITIONS
!
!  COMMON          NONE
!  BLOCKS
!
!  I/O             NONE
!
!  PRECISION       SINGLE
!
!  SPECIALIST      ROLAND SWEET
!
!  LANGUAGE        FORTRAN
!
!  HISTORY         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
!                  NATIONAL BUREAU OF STANDARDS (BOULDER).
!
!  ALGORITHM       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
!                  OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
!
!  PORTABILITY     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
!                  THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
!                  THE FUNCTION PIMACH.
!
!  REQUIRED        COS,SIN
!  RESIDENT
!  ROUTINES
!
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!            Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!            pp. 51-83.
!***ROUTINES CALLED  VRFTF1
!***END PROLOGUE  VRFFTF
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION       r(mdimr,n)  ,rt(mdimr,n)    ,wsave(n+15)
!***FIRST EXECUTABLE STATEMENT  VRFFTF
  IF (n == 1) RETURN
  CALL vrftf1 (m,n,r,rt,mdimr,wsave(1),wsave(n+1))
  RETURN
END SUBROUTINE vrfftf


SUBROUTINE vrffti (n,wsave) 1,1
!***BEGIN PROLOGUE  VRFFTI
!***DATE WRITTEN   860701   (YYMMDD)
!***REVISION DATE  900509   (YYMMDD)
!***CATEGORY NO.  J1A1
!***KEYWORDS  FAST FOURIER TRANSFORM, REAL PERIODIC TRANSFORM,
!          MULTIPLE SEQUENCES
!***AUTHOR  SWEET, R.A. (NIST) AND LINDGREN, L.L. (NIST)
!***PURPOSE  Initialization for VRFFTF and VRFFTB.
!***DESCRIPTION
!
!  Subroutine VRFFTI initializes the array WSAVE which is used in
!  both VRFFTF and VRFFTB.  The prime factorization of N together with
!  a tabulation of certain trigonometric functions are computed and
!  stored in the array WSAVE.
!
!  Input Parameter
!
!  N       the length of the sequence to be transformed.  There is no
!       restriction on N.
!
!  Output Parameter
!
!  WSAVE   a work array which must be dimensioned at least N+15.
!       The same work array can be used for both VRFFTF and VRFFTB
!       as long as N remains unchanged.  Different WSAVE arrays
!       are required for different values of N.  The contents of
!       WSAVE must not be changed between calls of VRFFTF or VRFFTB.
!
!
!           * * * * * * * * * * * * * * * * * * * * *
!           *                                       *
!           *         PROGRAM SPECIFICATIONS        *
!           *                                       *
!           * * * * * * * * * * * * * * * * * * * * *
!
!
!  DIMENSION OF    R(MDIMR,N), RT(MDIMR,N), WSAVE(N+15)
!  ARGUMENTS
!
!  LATEST          AUGUST 1, 1985
!  REVISION
!
!  SUBPROGRAMS     VRFFTI, VRFTI1, VRFFTF, VRFTF1, VRADF2, VRADF3,
!  REQUIRED        VRADF4, VRADF5, VRADFG, VRFFTB, VRFTB1, VRADB2,
!                  VRADB3, VRADB4, VRADB5, VRADBG, PIMACH
!
!  SPECIAL         NONE
!  CONDITIONS
!
!  COMMON          NONE
!  BLOCKS
!
!  I/O             NONE
!
!  PRECISION       SINGLE
!
!  SPECIALIST      ROLAND SWEET
!
!  LANGUAGE        FORTRAN
!
!  HISTORY         WRITTEN BY LINDA LINDGREN AND ROLAND SWEET AT THE
!                  NATIONAL BUREAU OF STANDARDS (BOULDER).
!
!  ALGORITHM       A REAL VARIANT OF THE STOCKHAM AUTOSORT VERSION
!                  OF THE COOLEY-TUKEY FAST FOURIER TRANSFORM.
!
!  PORTABILITY     AMERICAN NATIONAL STANDARDS INSTITUTE FORTRAN 77.
!                  THE ONLY MACHINE DEPENDENT CONSTANT IS LOCATED IN
!                  THE FUNCTION PIMACH.
!
!  REQUIRED        COS,SIN
!  RESIDENT
!  ROUTINES
!
!
!***REFERENCES  P. N. Swarztrauber, Vectorizing the FFTs, in Parallel
!            Computations, (G. Rodrigue, ed.), Academic Press, 1982,
!            pp. 51-83.
!***ROUTINES CALLED  VRFTI1
!***END PROLOGUE  VRFFTI
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION       wsave(n+15)
!***FIRST EXECUTABLE STATEMENT  VRFFTI
  IF (n == 1) RETURN
  CALL vrfti1 (n,wsave(1),wsave(n+1))
  RETURN
END SUBROUTINE vrffti


SUBROUTINE vrftb1 (m,n,c,ch,mdimc,wa,fac) 1,10
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION       ch(mdimc,n), c(mdimc,n), wa(n) ,fac(15)
  nf = fac(2)
  na = 0
  l1 = 1
  iw = 1
  DO k1=1,nf
    ip = fac(k1+2)
    l2 = ip*l1
    ido = n/l2
    idl1 = ido*l1
    IF (ip /= 4) GO TO 103
    ix2 = iw+ido
    ix3 = ix2+ido
    IF (na /= 0) GO TO 101
    CALL vradb4 (m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3))
    GO TO 102
    101    CALL vradb4 (m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3))
    102    na = 1-na
    GO TO 115
    103    IF (ip /= 2) GO TO 106
    IF (na /= 0) GO TO 104
    CALL vradb2 (m,ido,l1,c,ch,mdimc,wa(iw))
    GO TO 105
    104    CALL vradb2 (m,ido,l1,ch,c,mdimc,wa(iw))
    105    na = 1-na
    GO TO 115
    106    IF (ip /= 3) GO TO 109
    ix2 = iw+ido
    IF (na /= 0) GO TO 107
    CALL vradb3 (m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2))
    GO TO 108
    107    CALL vradb3 (m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2))
    108    na = 1-na
    GO TO 115
    109    IF (ip /= 5) GO TO 112
    ix2 = iw+ido
    ix3 = ix2+ido
    ix4 = ix3+ido
    IF (na /= 0) GO TO 110
    CALL vradb5 (m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    GO TO 111
    110 CALL vradb5 (m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    111    na = 1-na
    GO TO 115
    112    IF (na /= 0) GO TO 113
    CALL vradbg (m,ido,ip,l1,idl1,c,c,c,ch,ch,mdimc,wa(iw))
    GO TO 114
    113    CALL vradbg (m,ido,ip,l1,idl1,ch,ch,ch,c,c,mdimc,wa(iw))
    114    IF (ido == 1) na = 1-na
    115    l1 = l2
    iw = iw+(ip-1)*ido
  END DO
  scale=SQRT(1./n)
  IF (na == 0) GO TO 118
  DO j=1,n
    DO i=1,m
      c(i,j) = scale*ch(i,j)
    END DO
  END DO
  RETURN
  118 DO j=1,n
    DO i=1,m
      c(i,j)=scale*c(i,j)
    END DO
  END DO
  RETURN
END SUBROUTINE vrftb1


SUBROUTINE vrftf1 (m,n,c,ch,mdimc,wa,fac) 1,10
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION       ch(mdimc,n) ,c(mdimc,n)  ,wa(n)   ,fac(15)
  nf = fac(2)
  na = 1
  l2 = n
  iw = n
  DO k1=1,nf
    kh = nf-k1
    ip = fac(kh+3)
    l1 = l2/ip
    ido = n/l2
    idl1 = ido*l1
    iw = iw-(ip-1)*ido
    na = 1-na
    IF (ip /= 4) GO TO 102
    ix2 = iw+ido
    ix3 = ix2+ido
    IF (na /= 0) GO TO 101
    CALL vradf4 (m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3))
    GO TO 110
    101    CALL vradf4 (m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3))
    GO TO 110
    102    IF (ip /= 2) GO TO 104
    IF (na /= 0) GO TO 103
    CALL vradf2 (m,ido,l1,c,ch,mdimc,wa(iw))
    GO TO 110
    103    CALL vradf2 (m,ido,l1,ch,c,mdimc,wa(iw))
    GO TO 110
    104    IF (ip /= 3) GO TO 106
    ix2 = iw+ido
    IF (na /= 0) GO TO 105
    CALL vradf3 (m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2))
    GO TO 110
    105    CALL vradf3 (m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2))
    GO TO 110
    106    IF (ip /= 5) GO TO 108
    ix2 = iw+ido
    ix3 = ix2+ido
    ix4 = ix3+ido
    IF (na /= 0) GO TO 107
    CALL vradf5(m,ido,l1,c,ch,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    GO TO 110
    107 CALL vradf5 (m,ido,l1,ch,c,mdimc,wa(iw),wa(ix2),wa(ix3),wa(ix4))
    GO TO 110
    108    IF (ido == 1) na = 1-na
    IF (na /= 0) GO TO 109
    CALL vradfg (m,ido,ip,l1,idl1,c,c,c,ch,ch,mdimc,wa(iw))
    na = 1
    GO TO 110
    109    CALL vradfg (m,ido,ip,l1,idl1,ch,ch,ch,c,c,mdimc,wa(iw))
    na = 0
    110    l2 = l1
  END DO
  scale=SQRT(1./n)
  IF (na == 1) GO TO 113
  DO j=1,n
    DO i=1,m
      c(i,j) = scale*ch(i,j)
    END DO
  END DO
  RETURN
  113 DO j=1,n
    DO i=1,m
      c(i,j)=scale*c(i,j)
    END DO
  END DO
  RETURN
END SUBROUTINE vrftf1


SUBROUTINE vrfti1 (n,wa,fac) 1
!
!  VRFFTPK, VERSION 1, AUGUST 1985
!
  DIMENSION       wa(n)      ,fac(15)    ,ntryh(4)
  DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
  nl = n
  nf = 0
  j = 0
  101 j = j+1
  IF (j-4 > 0) THEN
    GO TO   103
  END IF
  ntry = ntryh(j)
  GO TO 104
  103 ntry = ntry+2
  104 nq = nl/ntry
  nr = nl-ntry*nq
  IF (nr == 0) THEN
    GO TO   105
  ELSE
    GO TO   101
  END IF
  105 nf = nf+1
  fac(nf+2) = ntry
  nl = nq
  IF (ntry /= 2) GO TO 107
  IF (nf == 1) GO TO 107
  DO i=2,nf
    ib = nf-i+2
    fac(ib+2) = fac(ib+1)
  END DO
  fac(3) = 2
  107 IF (nl /= 1) GO TO 104
  fac(1) = n
  fac(2) = nf
  tpi = 2.*pimach(1.0)
  argh = tpi/FLOAT(n)
  is = 0
  nfm1 = nf-1
  l1 = 1
  IF (nfm1 == 0) RETURN
  DO k1=1,nfm1
    ip = fac(k1+2)
    ld = 0
    l2 = l1*ip
    ido = n/l2
    ipm = ip-1
    DO j=1,ipm
      ld = ld+l1
      i = is
      argld = FLOAT(ld)*argh
      fi = 0.
      DO ii=3,ido,2
        i = i+2
        fi = fi+1.
        arg = fi*argld
        wa(i-1) = COS(arg)
        wa(i) = SIN(arg)
      END DO
      is = is+ido
    END DO
    l1 = l2
  END DO
  RETURN
END SUBROUTINE vrfti1