!                                               *****************
!                                               *   NFFT1.FOR   *
!                                               *  PURSER 1994  *
!                                               *****************
!
      BLOCKDATA DAT235
      COMMON/FFT235/ LN2,LN3,LN5
      DATA LN2/0/,LN3/0/,LN5/0/
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINES CFFT, DFFT                                     C
!                                                                              C
!   Fourier analyze (CFFT) or synthesize (DFFT) a line of complex data         C
!  RB     <-> Real part of data and transform                                  C
!  QB     <-> Imaginary part of data and transform                             C
!  N      --> Number of complex data of series (product of 2's,3's,5's)        C
!  PERIOD --> Period of one cycle                                              C
!  W      <-> Coefficient array of N reals (re)computed only when              C
!             JUMBLE(0)=0                                                      C
!  JUMBLE <-> Coefficient array of N integers (re)computed only when           C
!             JUMBLE(0) set to zero on input to this routine                   C
!                                                                              C
!   For related routines, see subroutine RUMBLE (in which the                  C
!   coefficient arrays are initialized), XCFFT, YCFFT (matrix versions         C
!   of CFFT), XDFFT, YDFFT (matrix versions of DFFT), RFFT (Fourier            C
!   analysis of real data), HFFT (Fourier synthesis of real data), and         C
!   matrix counterparts, XRFFT, YRFFT, XHFFT, YHFFT.                           C
!------------------------------------------------------------------------------C
      SUBROUTINE CFFT(RB,QB,N,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(0:*),QB(0:*),W(0:*),JUMBLE(0:*)
      RFAC=PERIOD/N
      GOTO 300
      ENTRY      DFFT(RB,QB,N,PERIOD,W,JUMBLE)
      RFAC=1./PERIOD

!  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
      DO J=N/2+1,N-1
       I=N-J
       T    =RB(J)
       RB(J)=RB(I)
       RB(I)=T
       T    =QB(J)
       QB(J)=QB(I)
       QB(I)=T
      ENDDO

300   IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

!  SCALE AND PERMUTE THE DATA:
      RB(0)=RB(0)*RFAC
      QB(0)=QB(0)*RFAC
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        T1=RB(I)
        RB(I)=RB(J)*RFAC
        RB(J)=T1
        T1=QB(I)
        QB(I)=QB(J)*RFAC
        QB(J)=T1
       ELSE
        RB(I)=RB(I)*RFAC
        QB(I)=QB(I)*RFAC
       ENDIF
      ENDDO

!  TRANSFORM THE DATA:
      MA=1
      MB=N
!  RADIX 4
      LS=0
      DO L=2,LN2,2
      MB=MB/4
      LS=L
      MA4=MA*4
      MB2=MB*2
      DO J=0,MA-1
      JMB=J*MB
      JMB2=J*MB2
      RF1=W(JMB)
      QF1=W(NH+JMB)
      RF2=W(JMB2)
      QF2=W(NH+JMB2)
      RF3=RF1*RF2-QF1*QF2
      QF3=RF1*QF2+QF1*RF2
      DO I=0,NM,MA4
      K0=I+J
      K1=K0+MA
      K2=K1+MA
      K3=K2+MA
      R0=RB(K0)
      R1=RB(K1)
      R2=RB(K2)
      R3=RB(K3)
      Q0=QB(K0)
      Q1=QB(K1)
      Q2=QB(K2)
      Q3=QB(K3)
      T1=R3*RF3-Q3*QF3  ! Q13
      Q3=R3*QF3+Q3*RF3  ! R13
      R3=R2*RF1-Q2*QF1  ! R12
      R2=R2*QF1+Q2*RF1  ! Q12
      Q2=Q3-R2          ! R23
      R2=Q3+R2          ! Q22
      Q3=R3+T1          ! R22
      T1=R3-T1          ! Q23
      R3=R1*RF2-Q1*QF2  ! R11
      Q1=R1*QF2+Q1*RF2  ! Q11
      R1=R0-R3          ! R21
      R0=R0+R3          ! R20
      RB(K3)=R1-Q2      ! R3
      RB(K1)=R1+Q2      ! R1
      Q2=Q0+Q1          ! Q20
      Q1=Q0-Q1          ! Q21
      QB(K0)=Q2+R2      ! Q0
      QB(K2)=Q2-R2      ! Q2
      RB(K2)=R0-Q3      ! R2
      RB(K0)=R0+Q3      ! R0
      QB(K3)=Q1-T1      ! Q3
      QB(K1)=Q1+T1      ! Q1
      ENDDO
      ENDDO
      MA=MA4
      ENDDO
      IF(LS.NE.LN2)THEN
!  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
      JMB=J*MB
      DO I=0,NM,MA2
      K0=J+I
      K1=K0+MA
      RF1=W(JMB)
      QF1=W(NH+JMB)
      R0=RB(K0)
      Q0=QB(K0)
      R1=RB(K1)
      Q1=QB(K1)
      T1=R1*QF1+Q1*RF1 ! Q11
      Q1=R1*RF1-Q1*QF1 ! R11
      RB(K1)=R0-Q1     ! R1
      RB(K0)=R0+Q1     ! R0
      QB(K1)=Q0-T1     ! Q1
      QB(K0)=Q0+T1     ! Q0
      ENDDO
      ENDDO
      MA=MA2
      ENDIF
!  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
      MB=MB/3
      MA3=MA*3
      DO J=0,MA-1
      JMB=J*MB
      RF1=W(JMB)
      QF1=W(NH+JMB)
      RF2=RF1*RF1-QF1*QF1
      QF2=2*RF1*QF1
      DO I=0,NM,MA3
      K0=I+J
      K1=K0+MA
      K2=K1+MA
      R1=RB(K1)
      Q1=QB(K1)
      R2=RB(K2)
      Q2=QB(K2)
      T1=R2*QF2+Q2*RF2  ! R12
      Q2=R2*RF2-Q2*QF2  ! Q12
      R2=R1*QF1+Q1*RF1  ! Q11
      R1=R1*RF1-Q1*QF1  ! R11
      Q1=R2+T1          ! Q21
      R2=(R2-T1)*QEP    ! R22
      T1=R1+Q2          ! R21
      R1=(R1-Q2)*QEP    ! Q22
      RB(K0)=RB(K0)+T1  ! R0
      QB(K0)=QB(K0)+Q1  ! Q0
      T1=RB(K0)-T1*REC  ! R21
      Q1=QB(K0)-Q1*REC  ! Q21
      QB(K2)=Q1-R1      ! Q2
      QB(K1)=Q1+R1      ! Q1
      RB(K1)=T1-R2      ! R1
      RB(K2)=T1+R2      ! R2
      ENDDO
      ENDDO
      MA=MA3
      ENDDO

      IF(LN5.GT.0)THEN
!  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
      MB=MB/5
      MA5=MA*5
      DO J=0,MA-1
      JMB=J*MB
      JMB2=JMB*2
      RF1=W(JMB)
      QF1=W(NH+JMB)
      RF2=W(JMB2)
      QF2=W(NH+JMB2)
      RF3=RF1*RF2-QF1*QF2
      QF3=RF1*QF2+QF1*RF2
      RF4=RF2*RF2-QF2*QF2
      QF4=2*RF2*QF2
      DO I=0,NM,MA5
      K0=I+J
      K1=K0+MA
      K2=K1+MA
      K3=K2+MA
      K4=K3+MA
      R1=RB(K1)
      R2=RB(K2)
      R3=RB(K3)
      R4=RB(K4)
      Q1=QB(K1)
      Q2=QB(K2)
      Q3=QB(K3)
      Q4=QB(K4)
      T1=R1*QF1+Q1*RF1        ! Q11
      R1=R1*RF1-Q1*QF1        ! R11
      Q1=R4*RF4-Q4*QF4        ! Q14
      R4=R4*QF4+Q4*RF4        ! R14
      Q4=R1-Q1                ! Q24
      R1=R1+Q1                ! R21
      Q1=T1+R4                ! Q21
      R4=T1-R4                ! R24
      T1=R3*RF3-Q3*QF3        ! Q13
      R3=R3*QF3+Q3*RF3        ! R13
      Q3=R2*QF2+Q2*RF2        ! Q12
      R2=R2*RF2-Q2*QF2        ! R12
      Q2=Q3+R3                ! Q22
      R3=Q3-R3                ! R23
      Q3=R2-T1                ! Q23
      R2=R2+T1                ! R22
      RB(K0)=RB(K0)+R1+R2     ! R0
      QB(K0)=QB(K0)+Q1+Q2     ! Q0
      T1=       R4*QZE+R3*QET ! R34
      R3=       R3*QZE-R4*QET ! R33
      R4=RB(K0)-R2*RZC-R1*REC ! R32
      R1=RB(K0)-R1*RZC-R2*REC ! R31
      RB(K2)=R4+R3            ! R2
      RB(K3)=R4-R3            ! R3
      RB(K4)=R1+T1            ! R4
      RB(K1)=R1-T1            ! R1
      T1=QB(K0)-Q1*RZC-Q2*REC ! Q31
      Q2=QB(K0)-Q2*RZC-Q1*REC ! Q32
      Q1=       Q3*QZE-Q4*QET ! Q33
      Q4=       Q4*QZE+Q3*QET ! Q34
      QB(K3)=Q2+Q1            ! Q3
      QB(K2)=Q2-Q1            ! Q2
      QB(K1)=T1+Q4            ! Q1
      QB(K4)=T1-Q4            ! Q4
      ENDDO
      ENDDO
      MA=MA5
      ENDDO
      ENDIF
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE GET235                                          C
!                                                                              C
!   Analyze the facorization of N in terms of 2's, 3's, and 5's and verify     C
!   that these are the only prime factors. Set LN2, LN3, LN5 respectively to   C
!   the number of powers of 2, 3, 5.                                           C
!  N      --> Number of data along the line of Fourier transformation          C
!------------------------------------------------------------------------------C
      SUBROUTINE GET235(N)
      COMMON/FFT235/ LN2,LN3,LN5
      NKT=N
      NM=N-1
      NH=N/2
      LN2=-1
400   NK=NKT
      NKT=NK/2
      LN2=LN2+1
      IF(2*NKT.EQ.NK)GOTO 400
      NKT=NK
      LN3=-1
401   NK=NKT
      NKT=NK/3
      LN3=LN3+1
      IF(3*NKT.EQ.NK)GOTO 401
      NKT=NK
      LN5=-1
402   NK=NKT
      NKT=NK/5
      LN5=LN5+1
      IF(5*NKT.EQ.NK)GOTO 402
      IF(NK.NE.1)STOP'N CANNOT BE FACTORED INTO ONLY 2s, 3s, 5s'
      RETURN
      END

!------------------------------------------------------------------------------C
!   R.J.Purser, National Meteorological Center, Washington D.C.  1994          C
!                   SUBROUTINE RUMBLE                                          C
!                                                                              C
!   Initialize coefficient arrays JUMBLE and TUMBLE for use with the           C
!   fast-Fourier-transform routines when the number N of data has the prime-   C
!   factorization 2**LN2*3**LN3*5**LN5                                         C
!                                                                              C
!  JUMBLE <-- Permutation of N data encoded as a sequence of transpositions    C
!  TUMBLE <-- Trigonometric coefficients for use in FFT. The first half are    C
!             the cosines, the second half the sines, of uniformly increasing  C
!             relevant angles                                                  C
!------------------------------------------------------------------------------C
      SUBROUTINE RUMBLE(JUMBLE,TUMBLE)
      PARAMETER(ML=20)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION ND(ML),MD(ML),JUMBLE(0:*),TUMBLE(0:*)
      LN=LN2+LN3+LN5
      N=2**LN2*3**LN3*5**LN5
      NM=N-1
      NH=N/2
      PI2ON=8.*ATAN(1.)/N
      DO I=0,NH-1
      ANG=PI2ON*I
      TUMBLE(I)=COS(ANG)
      TUMBLE(I+NH)=SIN(ANG)
      ENDDO
      ID=1
      IS=0
      DO I=1,LN5
       IS=IS+1
       MD(IS)=ID
       ID=ID*5
      ENDDO
      DO I=1,LN3
       IS=IS+1
       MD(IS)=ID
       ID=ID*3
      ENDDO
      DO I=1,LN2
       IS=IS+1
       MD(IS)=ID
       ID=ID*2
      ENDDO
      ID=1
      DO I=1,LN2
       ND(IS)=ID
       ID=ID*2
       IS=IS-1
      ENDDO
      DO I=1,LN3
       ND(IS)=ID
       ID=ID*3
       IS=IS-1
      ENDDO
      DO I=1,LN5
       ND(IS)=ID
       ID=ID*5
       IS=IS-1
      ENDDO
      JUMBLE(0)=N
      DO I=1,NM
       IR=I
       J=0
       DO L=1,LN
        KD=IR/ND(L)
        IR=IR-KD*ND(L)
        J=J+KD*MD(L)
       ENDDO
       JUMBLE(I)=J
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.LT.I)THEN
400        J=JUMBLE(J)
           IF(J.LT.I)GOTO 400
        JUMBLE(I)=J
       ENDIF
      ENDDO
      RETURN
      END

      SUBROUTINE XCFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(0:NCB-1,*),QB(0:NCB-1,*),W(0:*),JUMBLE(0:*)
      RFAC=PERIOD/N
!      GOTO 300
!
!      ENTRY     XDFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE)
!      RFAC=1./PERIOD
!
!C  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
!      DO J=N/2+1,N-1
!       I=N-J
!       DO IY=1,NY
!       T       =RB(J,IY)
!       RB(J,IY)=RB(I,IY)
!       RB(I,IY)=T
!       T       =QB(J,IY)
!       QB(J,IY)=QB(I,IY)
!       QB(I,IY)=T
!       ENDDO
!      ENDDO
!
300   IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

!  SCALE AND PERMUTE THE DATA:
      DO IY=1,NY
       RB(0,IY)=RB(0,IY)*RFAC
       QB(0,IY)=QB(0,IY)*RFAC
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        DO IY=1,NY
         T=RB(I,IY)
         RB(I,IY)=RB(J,IY)*RFAC
         RB(J,IY)=T
         T=QB(I,IY)
         QB(I,IY)=QB(J,IY)*RFAC
         QB(J,IY)=T
        ENDDO
       ELSE
        DO IY=1,NY
         RB(I,IY)=RB(I,IY)*RFAC
         QB(I,IY)=QB(I,IY)*RFAC
        ENDDO
       ENDIF
      ENDDO
!  TRANSFORM THE DATA:
      MA=1
      MB=N
!  RADIX 4
      LS=0
      DO L=2,LN2,2
       MB=MB/4
       LS=L
       MA4=MA*4
       MB2=MB*2
       DO J=0,MA-1
        JMB=J*MB
        JMB2=J*MB2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        DO I=0,NM,MA4
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         DO IY=1,NY
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3    ! Q13
          QB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3    ! R13
          RB(K3,IY)=RB(K2,IY)*RF1-QB(K2,IY)*QF1    ! R12
          RB(K2,IY)=RB(K2,IY)*QF1+QB(K2,IY)*RF1    ! Q12
          QB(K2,IY)=QB(K3,IY)-RB(K2,IY)            ! R23
          RB(K2,IY)=QB(K3,IY)+RB(K2,IY)            ! Q22
          QB(K3,IY)=RB(K3,IY)+T                    ! R22
          T        =RB(K3,IY)-T                    ! Q23
          RB(K3,IY)=RB(K1,IY)*RF2-QB(K1,IY)*QF2    ! R11
          QB(K1,IY)=RB(K1,IY)*QF2+QB(K1,IY)*RF2    ! Q11
          RB(K1,IY)=RB(K0,IY)-RB(K3,IY)            ! R21
          RB(K0,IY)=RB(K0,IY)+RB(K3,IY)            ! R20
          RB(K3,IY)=RB(K1,IY)-QB(K2,IY)            ! R3
          RB(K1,IY)=RB(K1,IY)+QB(K2,IY)            ! R1
          QB(K2,IY)=QB(K0,IY)+QB(K1,IY)            ! Q20
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)            ! Q21
          QB(K0,IY)=QB(K2,IY)+RB(K2,IY)            ! Q0
          QB(K2,IY)=QB(K2,IY)-RB(K2,IY)            ! Q2
          RB(K2,IY)=RB(K0,IY)-QB(K3,IY)            ! R2
          RB(K0,IY)=RB(K0,IY)+QB(K3,IY)            ! R0
          QB(K3,IY)=QB(K1,IY)-T                    ! Q3
          QB(K1,IY)=QB(K1,IY)+T                    ! Q1
         ENDDO
        ENDDO
       ENDDO
       MA=MA4
      ENDDO

      IF(LS.NE.LN2)THEN
!  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
       JMB=J*MB
       DO I=0,NM,MA2
        K0=J+I
        K1=K0+MA
        RF1=W(JMB)
        QF1=W(NH+JMB)
        DO IY=1,NY
         T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11
         QB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11
         RB(K1,IY)=RB(K0,IY)-QB(K1,IY)         ! R1
         RB(K0,IY)=RB(K0,IY)+QB(K1,IY)         ! R0
         QB(K1,IY)=QB(K0,IY)-T                 ! Q1
         QB(K0,IY)=QB(K0,IY)+T                 ! Q0
        ENDDO
       ENDDO
      ENDDO
      MA=MA2
      ENDIF
!  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
       MB=MB/3
       MA3=MA*3
       DO J=0,MA-1
        JMB=J*MB
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=RF1*RF1-QF1*QF1
        QF2=2*RF1*QF1
        DO I=0,NM,MA3
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         DO IY=1,NY
          T        =RB(K2,IY)*QF2+QB(K2,IY)*RF2     ! R12
          QB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2     ! Q12
          RB(K2,IY)=RB(K1,IY)*QF1+QB(K1,IY)*RF1     ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1     ! R11
          QB(K1,IY)=RB(K2,IY)+T                     ! Q21
          RB(K2,IY)=(RB(K2,IY)-T)*QEP               ! R22
          T        =RB(K1,IY)+QB(K2,IY)             ! R21
          RB(K1,IY)=(RB(K1,IY)-QB(K2,IY))*QEP       ! Q22
          RB(K0,IY)=RB(K0,IY)+T                     ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)             ! Q0
          T        =RB(K0,IY)-T*REC                 ! R21
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)*REC         ! Q21
          QB(K2,IY)=QB(K1,IY)-RB(K1,IY)             ! Q2
          QB(K1,IY)=QB(K1,IY)+RB(K1,IY)             ! Q1
          RB(K1,IY)=T        -RB(K2,IY)             ! R1
          RB(K2,IY)=T        +RB(K2,IY)             ! R2
         ENDDO
        ENDDO
       ENDDO
       MA=MA3
      ENDDO

      IF(LN5.GT.0)THEN
!  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
       MB=MB/5
       MA5=MA*5
       DO J=0,MA-1
        JMB=J*MB
        JMB2=JMB*2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        RF4=RF2*RF2-QF2*QF2
        QF4=2*RF2*QF2
        DO I=0,NM,MA5
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         K4=K3+MA
         DO IY=1,NY
          T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1           ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1           ! R11
          QB(K1,IY)=RB(K4,IY)*RF4-QB(K4,IY)*QF4           ! Q14
          RB(K4,IY)=RB(K4,IY)*QF4+QB(K4,IY)*RF4           ! R14
          QB(K4,IY)=RB(K1,IY)-QB(K1,IY)                   ! Q24
          RB(K1,IY)=RB(K1,IY)+QB(K1,IY)                   ! R21
          QB(K1,IY)=T        +RB(K4,IY)                   ! Q21
          RB(K4,IY)=T        -RB(K4,IY)                   ! R24
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3           ! Q13
          RB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3           ! R13
          QB(K3,IY)=RB(K2,IY)*QF2+QB(K2,IY)*RF2           ! Q12
          RB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2           ! R12
          QB(K2,IY)=QB(K3,IY)+RB(K3,IY)                   ! Q22
          RB(K3,IY)=QB(K3,IY)-RB(K3,IY)                   ! R23
          QB(K3,IY)=RB(K2,IY)-T                           ! Q23
          RB(K2,IY)=RB(K2,IY)+T                           ! R22
          RB(K0,IY)=RB(K0,IY)+RB(K1,IY)+RB(K2,IY)         ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)+QB(K2,IY)         ! Q0
          T        =RB(K4,IY)*QZE+RB(K3,IY)*QET           ! R34
          RB(K3,IY)=RB(K3,IY)*QZE-RB(K4,IY)*QET           ! R33
          RB(K4,IY)=RB(K0,IY)-RB(K2,IY)*RZC-RB(K1,IY)*REC ! R32
          RB(K1,IY)=RB(K0,IY)-RB(K1,IY)*RZC-RB(K2,IY)*REC ! R31
          RB(K2,IY)=RB(K4,IY)+RB(K3,IY)                   ! R2
          RB(K3,IY)=RB(K4,IY)-RB(K3,IY)                   ! R3
          RB(K4,IY)=RB(K1,IY)+T                           ! R4
          RB(K1,IY)=RB(K1,IY)-T                           ! R1
          T        =QB(K0,IY)-QB(K1,IY)*RZC-QB(K2,IY)*REC ! Q31
          QB(K2,IY)=QB(K0,IY)-QB(K2,IY)*RZC-QB(K1,IY)*REC ! Q32
          QB(K1,IY)=QB(K3,IY)*QZE-QB(K4,IY)*QET           ! Q33
          QB(K4,IY)=QB(K4,IY)*QZE+QB(K3,IY)*QET           ! Q34
          QB(K3,IY)=QB(K2,IY)+QB(K1,IY)                   ! Q3
          QB(K2,IY)=QB(K2,IY)-QB(K1,IY)                   ! Q2
          QB(K1,IY)=T        +QB(K4,IY)                   ! Q1
          QB(K4,IY)=T        -QB(K4,IY)                   ! Q4
         ENDDO
        ENDDO
       ENDDO
       MA=MA5
      ENDDO
      ENDIF
      RETURN
      END

      SUBROUTINE XDFFT(RB,QB,N,NY,NCB,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(0:NCB-1,*),QB(0:NCB-1,*),W(0:*),JUMBLE(0:*)
      RFAC=1./PERIOD

!  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
      DO J=N/2+1,N-1
       I=N-J
       DO IY=1,NY
        T       =RB(J,IY)
        RB(J,IY)=RB(I,IY)
        RB(I,IY)=T
        T       =QB(J,IY)
        QB(J,IY)=QB(I,IY)
        QB(I,IY)=T
       ENDDO
      ENDDO

300   IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

!  SCALE AND PERMUTE THE DATA:
      DO IY=1,NY
       RB(0,IY)=RB(0,IY)*RFAC
       QB(0,IY)=QB(0,IY)*RFAC
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        DO IY=1,NY
         T=RB(I,IY)
         RB(I,IY)=RB(J,IY)*RFAC
         RB(J,IY)=T
         T=QB(I,IY)
         QB(I,IY)=QB(J,IY)*RFAC
         QB(J,IY)=T
        ENDDO
       ELSE
        DO IY=1,NY
         RB(I,IY)=RB(I,IY)*RFAC
         QB(I,IY)=QB(I,IY)*RFAC
        ENDDO
       ENDIF
      ENDDO
!  TRANSFORM THE DATA:
      MA=1
      MB=N
!  RADIX 4
      LS=0
      DO L=2,LN2,2
       MB=MB/4
       LS=L
       MA4=MA*4
       MB2=MB*2
       DO J=0,MA-1
        JMB=J*MB
        JMB2=J*MB2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        DO I=0,NM,MA4
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         DO IY=1,NY
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3    ! Q13
          QB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3    ! R13
          RB(K3,IY)=RB(K2,IY)*RF1-QB(K2,IY)*QF1    ! R12
          RB(K2,IY)=RB(K2,IY)*QF1+QB(K2,IY)*RF1    ! Q12
          QB(K2,IY)=QB(K3,IY)-RB(K2,IY)            ! R23
          RB(K2,IY)=QB(K3,IY)+RB(K2,IY)            ! Q22
          QB(K3,IY)=RB(K3,IY)+T                    ! R22
          T        =RB(K3,IY)-T                    ! Q23
          RB(K3,IY)=RB(K1,IY)*RF2-QB(K1,IY)*QF2    ! R11
          QB(K1,IY)=RB(K1,IY)*QF2+QB(K1,IY)*RF2    ! Q11
          RB(K1,IY)=RB(K0,IY)-RB(K3,IY)            ! R21
          RB(K0,IY)=RB(K0,IY)+RB(K3,IY)            ! R20
          RB(K3,IY)=RB(K1,IY)-QB(K2,IY)            ! R3
          RB(K1,IY)=RB(K1,IY)+QB(K2,IY)            ! R1
          QB(K2,IY)=QB(K0,IY)+QB(K1,IY)            ! Q20
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)            ! Q21
          QB(K0,IY)=QB(K2,IY)+RB(K2,IY)            ! Q0
          QB(K2,IY)=QB(K2,IY)-RB(K2,IY)            ! Q2
          RB(K2,IY)=RB(K0,IY)-QB(K3,IY)            ! R2
          RB(K0,IY)=RB(K0,IY)+QB(K3,IY)            ! R0
          QB(K3,IY)=QB(K1,IY)-T                    ! Q3
          QB(K1,IY)=QB(K1,IY)+T                    ! Q1
         ENDDO
        ENDDO
       ENDDO
       MA=MA4
      ENDDO

      IF(LS.NE.LN2)THEN
!  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
       JMB=J*MB
       DO I=0,NM,MA2
        K0=J+I
        K1=K0+MA
        RF1=W(JMB)
        QF1=W(NH+JMB)
        DO IY=1,NY
         T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1 ! Q11
         QB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1 ! R11
         RB(K1,IY)=RB(K0,IY)-QB(K1,IY)         ! R1
         RB(K0,IY)=RB(K0,IY)+QB(K1,IY)         ! R0
         QB(K1,IY)=QB(K0,IY)-T                 ! Q1
         QB(K0,IY)=QB(K0,IY)+T                 ! Q0
        ENDDO
       ENDDO
      ENDDO
      MA=MA2
      ENDIF
!  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
       MB=MB/3
       MA3=MA*3
       DO J=0,MA-1
        JMB=J*MB
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=RF1*RF1-QF1*QF1
        QF2=2*RF1*QF1
        DO I=0,NM,MA3
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         DO IY=1,NY
          T        =RB(K2,IY)*QF2+QB(K2,IY)*RF2     ! R12
          QB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2     ! Q12
          RB(K2,IY)=RB(K1,IY)*QF1+QB(K1,IY)*RF1     ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1     ! R11
          QB(K1,IY)=RB(K2,IY)+T                     ! Q21
          RB(K2,IY)=(RB(K2,IY)-T)*QEP               ! R22
          T        =RB(K1,IY)+QB(K2,IY)             ! R21
          RB(K1,IY)=(RB(K1,IY)-QB(K2,IY))*QEP       ! Q22
          RB(K0,IY)=RB(K0,IY)+T                     ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)             ! Q0
          T        =RB(K0,IY)-T*REC                 ! R21
          QB(K1,IY)=QB(K0,IY)-QB(K1,IY)*REC         ! Q21
          QB(K2,IY)=QB(K1,IY)-RB(K1,IY)             ! Q2
          QB(K1,IY)=QB(K1,IY)+RB(K1,IY)             ! Q1
          RB(K1,IY)=T        -RB(K2,IY)             ! R1
          RB(K2,IY)=T        +RB(K2,IY)             ! R2
         ENDDO
        ENDDO
       ENDDO
       MA=MA3
      ENDDO

      IF(LN5.GT.0)THEN
!  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
       MB=MB/5
       MA5=MA*5
       DO J=0,MA-1
        JMB=J*MB
        JMB2=JMB*2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        RF4=RF2*RF2-QF2*QF2
        QF4=2*RF2*QF2
        DO I=0,NM,MA5
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         K4=K3+MA
         DO IY=1,NY
          T        =RB(K1,IY)*QF1+QB(K1,IY)*RF1           ! Q11
          RB(K1,IY)=RB(K1,IY)*RF1-QB(K1,IY)*QF1           ! R11
          QB(K1,IY)=RB(K4,IY)*RF4-QB(K4,IY)*QF4           ! Q14
          RB(K4,IY)=RB(K4,IY)*QF4+QB(K4,IY)*RF4           ! R14
          QB(K4,IY)=RB(K1,IY)-QB(K1,IY)                   ! Q24
          RB(K1,IY)=RB(K1,IY)+QB(K1,IY)                   ! R21
          QB(K1,IY)=T        +RB(K4,IY)                   ! Q21
          RB(K4,IY)=T        -RB(K4,IY)                   ! R24
          T        =RB(K3,IY)*RF3-QB(K3,IY)*QF3           ! Q13
          RB(K3,IY)=RB(K3,IY)*QF3+QB(K3,IY)*RF3           ! R13
          QB(K3,IY)=RB(K2,IY)*QF2+QB(K2,IY)*RF2           ! Q12
          RB(K2,IY)=RB(K2,IY)*RF2-QB(K2,IY)*QF2           ! R12
          QB(K2,IY)=QB(K3,IY)+RB(K3,IY)                   ! Q22
          RB(K3,IY)=QB(K3,IY)-RB(K3,IY)                   ! R23
          QB(K3,IY)=RB(K2,IY)-T                           ! Q23
          RB(K2,IY)=RB(K2,IY)+T                           ! R22
          RB(K0,IY)=RB(K0,IY)+RB(K1,IY)+RB(K2,IY)         ! R0
          QB(K0,IY)=QB(K0,IY)+QB(K1,IY)+QB(K2,IY)         ! Q0
          T        =RB(K4,IY)*QZE+RB(K3,IY)*QET           ! R34
          RB(K3,IY)=RB(K3,IY)*QZE-RB(K4,IY)*QET           ! R33
          RB(K4,IY)=RB(K0,IY)-RB(K2,IY)*RZC-RB(K1,IY)*REC ! R32
          RB(K1,IY)=RB(K0,IY)-RB(K1,IY)*RZC-RB(K2,IY)*REC ! R31
          RB(K2,IY)=RB(K4,IY)+RB(K3,IY)                   ! R2
          RB(K3,IY)=RB(K4,IY)-RB(K3,IY)                   ! R3
          RB(K4,IY)=RB(K1,IY)+T                           ! R4
          RB(K1,IY)=RB(K1,IY)-T                           ! R1
          T        =QB(K0,IY)-QB(K1,IY)*RZC-QB(K2,IY)*REC ! Q31
          QB(K2,IY)=QB(K0,IY)-QB(K2,IY)*RZC-QB(K1,IY)*REC ! Q32
          QB(K1,IY)=QB(K3,IY)*QZE-QB(K4,IY)*QET           ! Q33
          QB(K4,IY)=QB(K4,IY)*QZE+QB(K3,IY)*QET           ! Q34
          QB(K3,IY)=QB(K2,IY)+QB(K1,IY)                   ! Q3
          QB(K2,IY)=QB(K2,IY)-QB(K1,IY)                   ! Q2
          QB(K1,IY)=T        +QB(K4,IY)                   ! Q1
          QB(K4,IY)=T        -QB(K4,IY)                   ! Q4
         ENDDO
        ENDDO
       ENDDO
       MA=MA5
      ENDDO
      ENDIF
      RETURN
      END

      SUBROUTINE YCFFT(RB,QB,N,NX,NDX,NCB,PERIOD,W,JUMBLE)
      COMMON/FFT235/ LN2,LN3,LN5
      DIMENSION RB(NCB,0:*),QB(NCB,0:*),W(0:*),JUMBLE(0:*)
      RFAC=PERIOD/N
      GOTO 300

      ENTRY      YDFFT(RB,QB,N,NX,NDX,NCB,PERIOD,W,JUMBLE)
      RFAC=1./PERIOD

!  FOR FOURIER SYNTHESIS, REVERSE THE ORDER OF WAVENUMBERS:
      DO J=N/2+1,N-1
       I=N-J
       DO IX=1,NX,NDX
        T       =RB(IX,J)
        RB(IX,J)=RB(IX,I)
        RB(IX,I)=T
        T       =QB(IX,J)
        QB(IX,J)=QB(IX,I)
        QB(IX,I)=T
       ENDDO
      ENDDO

300   IF(N.NE.2**LN2*3**LN3*5**LN5)CALL GET235(N)
      IF(N.NE.JUMBLE(0))CALL RUMBLE(JUMBLE,W)
      NM=N-1
      NH=N/2

!  SCALE AND PERMUTE THE DATA:
      DO IX=1,NX,NDX
       RB(IX,0)=RB(IX,0)*RFAC
       QB(IX,0)=QB(IX,0)*RFAC
      ENDDO
      DO I=1,NM
       J=JUMBLE(I)
       IF(J.GT.I)THEN
        DO IX=1,NX,NDX
         T=RB(IX,I)
         RB(IX,I)=RB(IX,J)*RFAC
         RB(IX,J)=T
         T=QB(IX,I)
         QB(IX,I)=QB(IX,J)*RFAC
         QB(IX,J)=T
        ENDDO
       ELSE
        DO IX=1,NX,NDX
         RB(IX,I)=RB(IX,I)*RFAC
         QB(IX,I)=QB(IX,I)*RFAC
        ENDDO
       ENDIF
      ENDDO

!  TRANSFORM THE DATA:
      MA=1
      MB=N
!  RADIX 4
      LS=0
      DO L=2,LN2,2
       MB=MB/4
       LS=L
       MA4=MA*4
       MB2=MB*2
       DO J=0,MA-1
        JMB=J*MB
        JMB2=J*MB2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        DO I=0,NM,MA4
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         DO IX=1,NX,NDX
          T        =RB(IX,K3)*RF3-QB(IX,K3)*QF3     ! Q13
          QB(IX,K3)=RB(IX,K3)*QF3+QB(IX,K3)*RF3     ! R13
          RB(IX,K3)=RB(IX,K2)*RF1-QB(IX,K2)*QF1     ! R12
          RB(IX,K2)=RB(IX,K2)*QF1+QB(IX,K2)*RF1     ! Q12
          QB(IX,K2)=QB(IX,K3)-RB(IX,K2)             ! R23
          RB(IX,K2)=QB(IX,K3)+RB(IX,K2)             ! Q22
          QB(IX,K3)=RB(IX,K3)+T                     ! R22
          T        =RB(IX,K3)-T                     ! Q23
          RB(IX,K3)=RB(IX,K1)*RF2-QB(IX,K1)*QF2     ! R11
          QB(IX,K1)=RB(IX,K1)*QF2+QB(IX,K1)*RF2     ! Q11
          RB(IX,K1)=RB(IX,K0)-RB(IX,K3)             ! R21
          RB(IX,K0)=RB(IX,K0)+RB(IX,K3)             ! R20
          RB(IX,K3)=RB(IX,K1)-QB(IX,K2)             ! R3
          RB(IX,K1)=RB(IX,K1)+QB(IX,K2)             ! R1
          QB(IX,K2)=QB(IX,K0)+QB(IX,K1)             ! Q20
          QB(IX,K1)=QB(IX,K0)-QB(IX,K1)             ! Q21
          QB(IX,K0)=QB(IX,K2)+RB(IX,K2)             ! Q0
          QB(IX,K2)=QB(IX,K2)-RB(IX,K2)             ! Q2
          RB(IX,K2)=RB(IX,K0)-QB(IX,K3)             ! R2
          RB(IX,K0)=RB(IX,K0)+QB(IX,K3)             ! R0
          QB(IX,K3)=QB(IX,K1)-T                     ! Q3
          QB(IX,K1)=QB(IX,K1)+T                     ! Q1
         ENDDO
        ENDDO
       ENDDO
       MA=MA4
      ENDDO
      IF(LS.NE.LN2)THEN

!  RADIX 2
      MB=MB/2
      MA2=MA*2
      DO J=0,MA-1
       JMB=J*MB
       DO I=0,NM,MA2
        K0=J+I
        K1=K0+MA
        RF1=W(JMB)
        QF1=W(NH+JMB)
        DO IX=1,NX,NDX
         T        =RB(IX,K1)*QF1+QB(IX,K1)*RF1 ! Q11
         QB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1 ! R11
         RB(IX,K1)=RB(IX,K0)-QB(IX,K1)         ! R1
         RB(IX,K0)=RB(IX,K0)+QB(IX,K1)         ! R0
         QB(IX,K1)=QB(IX,K0)-T                 ! Q1
         QB(IX,K0)=QB(IX,K0)+T                 ! Q0
        ENDDO
       ENDDO
      ENDDO
      MA=MA2
      ENDIF

!  RADIX 3
      REP=-.5
      REC=1.5
      QEP=.5*SQRT(3.)
      DO L=1,LN3
       MB=MB/3
       MA3=MA*3
       DO J=0,MA-1
        JMB=J*MB
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=RF1*RF1-QF1*QF1
        QF2=2*RF1*QF1
        DO I=0,NM,MA3
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         T        =RB(IX,K2)*QF2+QB(IX,K2)*RF2     ! R12
         QB(IX,K2)=RB(IX,K2)*RF2-QB(IX,K2)*QF2     ! Q12
         RB(IX,K2)=RB(IX,K1)*QF1+QB(IX,K1)*RF1     ! Q11
         RB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1     ! R11
         QB(IX,K1)=RB(IX,K2)+T                     ! Q21
         RB(IX,K2)=(RB(IX,K2)-T        )*QEP       ! R22
         T        =RB(IX,K1)+QB(IX,K2)             ! R21
         RB(IX,K1)=(RB(IX,K1)-QB(IX,K2))*QEP       ! Q22
         RB(IX,K0)=RB(IX,K0)+T                     ! R0
         QB(IX,K0)=QB(IX,K0)+QB(IX,K1)             ! Q0
         T        =RB(IX,K0)-T        *REC         ! R21
         QB(IX,K1)=QB(IX,K0)-QB(IX,K1)*REC         ! Q21
         QB(IX,K2)=QB(IX,K1)-RB(IX,K1)             ! Q2
         QB(IX,K1)=QB(IX,K1)+RB(IX,K1)             ! Q1
         RB(IX,K1)=T        -RB(IX,K2)             ! R1
         RB(IX,K2)=T        +RB(IX,K2)             ! R2
        ENDDO
       ENDDO
       MA=MA3
      ENDDO
      IF(LN5.GT.0)THEN

!  RADIX 5
      NZE=N/5
      RZE=W(NZE)
      QZE=W(NH+NZE)
      RZC=1.-RZE
      RET=RZE*RZE-QZE*QZE
      QET=2*RZE*QZE
      REC=1.-RET
      DO L=1,LN5
       MB=MB/5
       MA5=MA*5
       DO J=0,MA-1
        JMB=J*MB
        JMB2=JMB*2
        RF1=W(JMB)
        QF1=W(NH+JMB)
        RF2=W(JMB2)
        QF2=W(NH+JMB2)
        RF3=RF1*RF2-QF1*QF2
        QF3=RF1*QF2+QF1*RF2
        RF4=RF2*RF2-QF2*QF2
        QF4=2*RF2*QF2
        DO I=0,NM,MA5
         K0=I+J
         K1=K0+MA
         K2=K1+MA
         K3=K2+MA
         K4=K3+MA
         T        =RB(IX,K1)*QF1+QB(IX,K1)*RF1           ! Q11
         RB(IX,K1)=RB(IX,K1)*RF1-QB(IX,K1)*QF1           ! R11
         QB(IX,K1)=RB(IX,K4)*RF4-QB(IX,K4)*QF4           ! Q14
         RB(IX,K4)=RB(IX,K4)*QF4+QB(IX,K4)*RF4           ! R14
         QB(IX,K4)=RB(IX,K1)-QB(IX,K1)                   ! Q24
         RB(IX,K1)=RB(IX,K1)+QB(IX,K1)                   ! R21
         QB(IX,K1)=T        +RB(IX,K4)                   ! Q21
         RB(IX,K4)=T        -RB(IX,K4)                   ! R24
         T        =RB(IX,K3)*RF3-QB(IX,K3)*QF3           ! Q13
         RB(IX,K3)=RB(IX,K3)*QF3+QB(IX,K3)*RF3           ! R13
         QB(IX,K3)=RB(IX,K2)*QF2+QB(IX,K2)*RF2           ! Q12
         RB(IX,K2)=RB(IX,K2)*RF2-QB(IX,K2)*QF2           ! R12
         QB(IX,K2)=QB(IX,K3)+RB(IX,K3)                   ! Q22
         RB(IX,K3)=QB(IX,K3)-RB(IX,K3)                   ! R23
         QB(IX,K3)=RB(IX,K2)-T                           ! Q23
         RB(IX,K2)=RB(IX,K2)+T                           ! R22
         RB(IX,K0)=RB(IX,K0)+RB(IX,K1)+RB(IX,K2)         ! R0
         QB(IX,K0)=QB(IX,K0)+QB(IX,K1)+QB(IX,K2)         ! Q0
         T        =RB(IX,K4)*QZE+RB(IX,K3)*QET           ! R34
         RB(IX,K3)=RB(IX,K3)*QZE-RB(IX,K4)*QET           ! R33
         RB(IX,K4)=RB(IX,K0)-RB(IX,K2)*RZC-RB(IX,K1)*REC ! R32
         RB(IX,K1)=RB(IX,K0)-RB(IX,K1)*RZC-RB(IX,K2)*REC ! R31
         RB(IX,K2)=RB(IX,K4)+RB(IX,K3)                   ! R2
         RB(IX,K3)=RB(IX,K4)-RB(IX,K3)                   ! R3
         RB(IX,K4)=RB(IX,K1)+T                           ! R4
         RB(IX,K1)=RB(IX,K1)-T                           ! R1
         T        =QB(IX,K0)-QB(IX,K1)*RZC-QB(IX,K2)*REC ! Q31
         QB(IX,K2)=QB(IX,K0)-QB(IX,K2)*RZC-QB(IX,K1)*REC ! Q32
         QB(IX,K1)=QB(IX,K3)*QZE-QB(IX,K4)*QET           ! Q33
         QB(IX,K4)=QB(IX,K4)*QZE+QB(IX,K3)*QET           ! Q34
         QB(IX,K3)=QB(IX,K2)+QB(IX,K1)                   ! Q3
         QB(IX,K2)=QB(IX,K2)-QB(IX,K1)                   ! Q2
         QB(IX,K1)=T        +QB(IX,K4)                   ! Q1
         QB(IX,K4)=T        -QB(IX,K4)                   ! Q4
        ENDDO
       ENDDO
       MA=MA5
      ENDDO
      ENDIF
      RETURN
      END
