C
C     ..................................................................
C
C        SUBROUTINE HARM
C
C        PURPOSE
C           PERFORMS DISCRETE COMPLEX FOURIER TRANSFORMS ON A COMPLEX
C           THREE DIMENSIONAL ARRAY
C
C        USAGE
C           CALL HARM (A,M,INV,S,IFSET,IFERR)
C
C        DESCRIPTION OF PARAMETERS
C           A     - AS INPUT, A CONTAINS THE COMPLEX, 3-DIMENSIONAL
C                   ARRAY TO BE TRANSFORMED.  THE REAL PART OF
C                   A(I1,I2,I3) IS STORED IN VECTOR FASHION IN A CELL
C                   WITH INDEX 2*(I3*N1*N2 + I2*N1 + I1) + 1 WHERE
C                   NI = 2**M(I), I=1,2,3 AND I1 = 0,1,...,N1-1 ETC.
C                   THE IMAGINARY PART IS IN THE CELL IMMEDIATELY
C                   FOLLOWING.  NOTE THAT THE SUBSCRIPT I1 INCREASES
C                   MOST RAPIDLY AND I3 INCREASES LEAST RAPIDLY.
C                   AS OUTPUT, A CONTAINS THE COMPLEX FOURIER
C                   TRANSFORM.  THE NUMBER OF CORE LOCATIONS OF
C                   ARRAY A IS 2*(N1*N2*N3)
C           M     - A THREE CELL VECTOR WHICH DETERMINES THE SIZES
C                   OF THE 3 DIMENSIONS OF THE ARRAY A.   THE SIZE,
C                   NI, OF THE I DIMENSION OF A IS 2**M(I), I = 1,2,3
C           INV   - A VECTOR WORK AREA FOR BIT AND INDEX MANIPULATION
C                   OF DIMENSION ONE FOURTH OF THE QUANTITY
C                   MAX(N1,N2,N3)
C           S     - A VECTOR WORK AREA FOR SINE TABLES WITH DIMENSION
C                   THE SAME AS INV
C           IFSET - AN OPTION PARAMETER WITH THE FOLLOWING SETTINGS
C                      0    SET UP SINE AND INV TABLES ONLY
C                      1    SET UP SINE AND INV TABLES ONLY AND
C                           CALCULATE FOURIER TRANSFORM
C                     -1    SET UP SINE AND INV TABLES ONLY AND
C                           CALCULATE INVERSE FOURIER TRANSFORM (FOR
C                           THE MEANING OF INVERSE SEE THE EQUATIONS
C                           UNDER METHOD BELOW)
C                      2    CALCULATE FOURIER TRANSFORM ONLY (ASSUME
C                           SINE AND INV TABLES EXIST)
C                     -2    CALCULATE INVERSE FOURIER TRANSFORM ONLY
C                           (ASSUME SINE AND INV TABLES EXIST)
C           IFERR - ERROR INDICATOR.   WHEN IFSET IS 0,+1,-1,
C                   IFERR = 1 MEANS THE MAXIMUM M(I) IS GREATER THAN
C                  20 , I=1,2,3   WHEN IFSET IS 2,-2 , IFERR = 1
C                   MEANS THAT THE SINE AND INV TABLES ARE NOT LARGE
C                   ENOUGH OR HAVE NOT BEEN COMPUTED .
C                   IF ON RETURN IFERR = 0 THEN NONE OF THE ABOVE
C                   CONDITIONS ARE PRESENT
C
C        REMARKS
C           THIS SUBROUTINE IS TO BE USED FOR COMPLEX, 3-DIMENSIONAL
C           ARRAYS IN WHICH EACH DIMENSION IS A POWER OF 2.  THE
C           MAXIMUM M(I) MUST NOT BE LESS THAN 3 OR GREATER THAN 20,
C           I = 1,2,3
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           FOR IFSET = +1, OR +2, THE FOURIER TRANSFORM OF COMPLEX
C           ARRAY A IS OBTAINED.
C
C                  N1-1   N2-1   N3-1                L1   L2   L3
C     X(J1,J2,J3)=SUM    SUM    SUM    A(K1,K2,K3)*W1  *W2  *W3
C                  K1=0   K2=0   K3=0
C
C                  WHERE WI IS THE N(I) ROOT OF UNITY AND L1=K1*J1,
C                        L2=K2*J2, L3=K3*J3
C
C
C           FOR IFSET = -1, OR -2, THE INVERSE FOURIER TRANSFORM A OF
C           COMPLEX ARRAY X IS OBTAINED.
C
C     A(K1,K2,K3)=
C               1      N1-1   N2-1   N3-1                -L1  -L2  -L3
C           -------- *SUM    SUM    SUM    X(J1,J2,J3)*W1  *W2  *W3
C           N1*N2*N3   J1=0   J2=0   J3=0
C
C
C           SEE J.W. COOLEY AND J.W. TUKEY, 'AN ALGORITHM FOR THE
C           MACHINE CALCULATION OF COMPLEX FOURIER SERIES',
C           MATHEMATICS OF COMPUTATIONS, VOL. 19 (APR. 1965), P. 297.
C
C     ..................................................................
C
      SUBROUTINE HARM(A,M,INV,S,IFSET, IFERR)
      DIMENSION A(1),INV(1),S(1),N(3),M(3),NP(3),W(2),W2(2),W3(2)
      EQUIVALENCE (N(1),N1),(N(2),N2),(N(3),N3)
   10 IF( IABS(IFSET) - 1) 900,900,12
   12 MTT=MAX0(M(1),M(2),M(3)) -2
      ROOT2 = SQRT(2.)
      IF (MTT-MT ) 14,14,13
   13 IFERR=1
      RETURN
   14 IFERR=0
      M1=M(1)
      M2=M(2)
      M3=M(3)
      N1=2**M1
      N2=2**M2
      N3=2**M3
   16 IF(IFSET) 18,18,20
   18 NX= N1*N2*N3
      FN = NX
      DO 19 I = 1,NX
      A(2*I-1) = A(2*I-1)/FN
   19 A(2*I) = -A(2*I)/FN
   20 NP(1)=N1*2
      NP(2)= NP(1)*N2
      NP(3)=NP(2)*N3
      DO 250 ID=1,3
      IL = NP(3)-NP(ID)
      IL1 = IL+1
      MI = M(ID)
      IF (MI)250,250,30
   30 IDIF=NP(ID)
      KBIT=NP(ID)
      MEV = 2*(MI/2)
      IF (MI - MEV )60,60,40
C
C     M IS ODD. DO L=1 CASE
   40 KBIT=KBIT/2
      KL=KBIT-2
      DO 50 I=1,IL1,IDIF
      KLAST=KL+I
      DO 50 K=I,KLAST,2
      KD=K+KBIT
C
C     DO ONE STEP WITH L=1,J=0
C     A(K)=A(K)+A(KD)
C     A(KD)=A(K)-A(KD)
C
      T=A(KD)
      A(KD)=A(K)-T
      A(K)=A(K)+T
      T=A(KD+1)
      A(KD+1)=A(K+1)-T
   50 A(K+1)=A(K+1)+T
      IF (MI - 1)250,250,52
   52 LFIRST =3
C
C     DEF - JLAST = 2**(L-2) -1
      JLAST=1
      GO TO 70
C
C     M IS EVEN
   60 LFIRST = 2
      JLAST=0
   70 DO 240 L=LFIRST,MI,2
      JJDIF=KBIT
      KBIT=KBIT/4
      KL=KBIT-2
C
C     DO FOR J=0
      DO 80 I=1,IL1,IDIF
      KLAST=I+KL
      DO 80 K=I,KLAST,2
      K1=K+KBIT
      K2=K1+KBIT
      K3=K2+KBIT
C
C     DO TWO STEPS WITH J=0
C     A(K)=A(K)+A(K2)
C     A(K2)=A(K)-A(K2)
C     A(K1)=A(K1)+A(K3)
C     A(K3)=A(K1)-A(K3)
C
C     A(K)=A(K)+A(K1)
C     A(K1)=A(K)-A(K1)
C     A(K2)=A(K2)+A(K3)*I
C     A(K3)=A(K2)-A(K3)*I
C
      T=A(K2)
      A(K2)=A(K)-T
      A(K)=A(K)+T
      T=A(K2+1)
      A(K2+1)=A(K+1)-T
      A(K+1)=A(K+1)+T
C
      T=A(K3)
      A(K3)=A(K1)-T
      A(K1)=A(K1)+T
      T=A(K3+1)
      A(K3+1)=A(K1+1)-T
      A(K1+1)=A(K1+1)+T
C
      T=A(K1)
      A(K1)=A(K)-T
      A(K)=A(K)+T
      T=A(K1+1)
      A(K1+1)=A(K+1)-T
      A(K+1)=A(K+1)+T
C
      R=-A(K3+1)
      T = A(K3)
      A(K3)=A(K2)-R
      A(K2)=A(K2)+R
      A(K3+1)=A(K2+1)-T
   80 A(K2+1)=A(K2+1)+T
      IF (JLAST) 235,235,82
   82 JJ=JJDIF   +1
C
C     DO FOR J=1
      ILAST= IL +JJ
      DO 85 I = JJ,ILAST,IDIF
      KLAST = KL+I
      DO 85 K=I,KLAST,2
      K1 = K+KBIT
      K2 = K1+KBIT
      K3 = K2+KBIT
C
C     LETTING W=(1+I)/ROOT2,W3=(-1+I)/ROOT2,W2=I,
C     A(K)=A(K)+A(K2)*I
C     A(K2)=A(K)-A(K2)*I
C     A(K1)=A(K1)*W+A(K3)*W3
C     A(K3)=A(K1)*W-A(K3)*W3
C
C     A(K)=A(K)+A(K1)
C     A(K1)=A(K)-A(K1)
C     A(K2)=A(K2)+A(K3)*I
C     A(K3)=A(K2)-A(K3)*I
C
      R =-A(K2+1)
      T = A(K2)
      A(K2) = A(K)-R
      A(K) = A(K)+R
      A(K2+1)=A(K+1)-T
      A(K+1)=A(K+1)+T
C
      AWR=A(K1)-A(K1+1)
      AWI = A(K1+1)+A(K1)
      R=-A(K3)-A(K3+1)
      T=A(K3)-A(K3+1)
      A(K3)=(AWR-R)/ROOT2
      A(K3+1)=(AWI-T)/ROOT2
      A(K1)=(AWR+R)/ROOT2
      A(K1+1)=(AWI+T)/ROOT2
      T= A(K1)
      A(K1)=A(K)-T
      A(K)=A(K)+T
      T=A(K1+1)
      A(K1+1)=A(K+1)-T
      A(K+1)=A(K+1)+T
      R=-A(K3+1)
      T=A(K3)
      A(K3)=A(K2)-R
      A(K2)=A(K2)+R
      A(K3+1)=A(K2+1)-T
   85 A(K2+1)=A(K2+1)+T
      IF(JLAST-1) 235,235,90
   90 JJ= JJ + JJDIF
C
C     NOW DO THE REMAINING J'S
      DO 230 J=2,JLAST
C
C     FETCH W'S
C     DEF- W=W**INV(J), W2=W**2, W3=W**3
   96 I=INV(J+1)
   98 IC=NT-I
      W(1)=S(IC)
      W(2)=S(I)
      I2=2*I
      I2C=NT-I2
      IF(I2C)120,110,100
C
C     2*I IS IN FIRST QUADRANT
  100 W2(1)=S(I2C)
      W2(2)=S(I2)
      GO TO 130
  110 W2(1)=0.
      W2(2)=1.
      GO TO 130
C
C     2*I IS IN SECOND QUADRANT
  120 I2CC = I2C+NT
      I2C=-I2C
      W2(1)=-S(I2C)
      W2(2)=S(I2CC)
  130 I3=I+I2
      I3C=NT-I3
      IF(I3C)160,150,140
C
C     I3 IN FIRST QUADRANT
  140 W3(1)=S(I3C)
      W3(2)=S(I3)
      GO TO 200
  150 W3(1)=0.
      W3(2)=1.
      GO TO 200
C
  160 I3CC=I3C+NT
      IF(I3CC)190,180,170
C
C     I3 IN SECOND QUADRANT
  170 I3C=-I3C
      W3(1)=-S(I3C)
      W3(2)=S(I3CC)
      GO TO 200
  180 W3(1)=-1.
      W3(2)=0.
      GO TO 200
C
C     3*I IN THIRD QUADRANT
  190 I3CCC=NT+I3CC
      I3CC = -I3CC
      W3(1)=-S(I3CCC)
      W3(2)=-S(I3CC)
  200 ILAST=IL+JJ
      DO 220 I=JJ,ILAST,IDIF
      KLAST=KL+I
      DO 220 K=I,KLAST,2
      K1=K+KBIT
      K2=K1+KBIT
      K3=K2+KBIT
C
C     DO TWO STEPS WITH J NOT 0
C     A(K)=A(K)+A(K2)*W2
C     A(K2)=A(K)-A(K2)*W2
C     A(K1)=A(K1)*W+A(K3)*W3
C     A(K3)=A(K1)*W-A(K3)*W3
C
C     A(K)=A(K)+A(K1)
C     A(K1)=A(K)-A(K1)
C     A(K2)=A(K2)+A(K3)*I
C     A(K3)=A(K2)-A(K3)*I
C
      R=A(K2)*W2(1)-A(K2+1)*W2(2)
      T=A(K2)*W2(2)+A(K2+1)*W2(1)
      A(K2)=A(K)-R
      A(K)=A(K)+R
      A(K2+1)=A(K+1)-T
      A(K+1)=A(K+1)+T
C
      R=A(K3)*W3(1)-A(K3+1)*W3(2)
      T=A(K3)*W3(2)+A(K3+1)*W3(1)
      AWR=A(K1)*W(1)-A(K1+1)*W(2)
      AWI=A(K1)*W(2)+A(K1+1)*W(1)
      A(K3)=AWR-R
      A(K3+1)=AWI-T
      A(K1)=AWR+R
      A(K1+1)=AWI+T
      T=A(K1)
      A(K1)=A(K)-T
      A(K)=A(K)+T
      T=A(K1+1)
      A(K1+1)=A(K+1)-T
      A(K+1)=A(K+1)+T
      R=-A(K3+1)
      T=A(K3)
      A(K3)=A(K2)-R
      A(K2)=A(K2)+R
      A(K3+1)=A(K2+1)-T
  220 A(K2+1)=A(K2+1)+T
C     END OF I AND K LOOPS
C
  230 JJ=JJDIF+JJ
C     END OF J-LOOP
C
  235 JLAST=4*JLAST+3
  240 CONTINUE
C     END OF  L  LOOP
C
  250 CONTINUE
C     END OF  ID  LOOP
C
C     WE NOW HAVE THE COMPLEX FOURIER SUMS BUT THEIR ADDRESSES ARE
C     BIT-REVERSED.  THE FOLLOWING ROUTINE PUTS THEM IN ORDER
      NTSQ=NT*NT
      M3MT=M3-MT
  350 IF(M3MT) 370,360,360
C
C     M3 GR. OR EQ. MT
  360 IGO3=1
      N3VNT=N3/NT
      MINN3=NT
      GO TO 380
C
C     M3 LESS THAN MT
  370 IGO3=2
      N3VNT=1
      NTVN3=NT/N3
      MINN3=N3
  380 JJD3 = NTSQ/N3
      M2MT=M2-MT
  450 IF (M2MT)470,460,460
C
C     M2 GR. OR EQ. MT
  460 IGO2=1
      N2VNT=N2/NT
      MINN2=NT
      GO TO 480
C
C     M2 LESS THAN MT
  470 IGO2 = 2
      N2VNT=1
      NTVN2=NT/N2
      MINN2=N2
  480 JJD2=NTSQ/N2
      M1MT=M1-MT
  550 IF(M1MT)570,560,560
C
C     M1 GR. OR EQ. MT
  560 IGO1=1
      N1VNT=N1/NT
      MINN1=NT
      GO TO 580
C
C     M1 LESS THAN MT
  570 IGO1=2
      N1VNT=1
      NTVN1=NT/N1
      MINN1=N1
  580 JJD1=NTSQ/N1
  600 JJ3=1
      J=1
      DO 880 JPP3=1,N3VNT
      IPP3=INV(JJ3)
      DO 870 JP3=1,MINN3
      GO TO (610,620),IGO3
  610 IP3=INV(JP3)*N3VNT
      GO TO 630
  620 IP3=INV(JP3)/NTVN3
  630 I3=(IPP3+IP3)*N2
  700 JJ2=1
      DO 870 JPP2=1,N2VNT
      IPP2=INV(JJ2)+I3
      DO 860 JP2=1,MINN2
      GO TO (710,720),IGO2
  710 IP2=INV(JP2)*N2VNT
      GO TO 730
  720 IP2=INV(JP2)/NTVN2
  730 I2=(IPP2+IP2)*N1
  800 JJ1=1
      DO 860 JPP1=1,N1VNT
      IPP1=INV(JJ1)+I2
      DO 850 JP1=1,MINN1
      GO TO (810,820),IGO1
  810 IP1=INV(JP1)*N1VNT
      GO TO 830
  820 IP1=INV(JP1)/NTVN1
  830 I=2*(IPP1+IP1)+1
      IF (J-I) 840,850,850
  840 T=A(I)
      A(I)=A(J)
      A(J)=T
      T=A(I+1)
      A(I+1)=A(J+1)
      A(J+1)=T
  850 J=J+2
  860 JJ1=JJ1+JJD1
C     END OF JPP1 AND JP2
C
  870 JJ2=JJ2+JJD2
C     END OF JPP2 AND JP3 LOOPS
C
  880 JJ3 = JJ3+JJD3
C     END OF JPP3 LOOP
C
  890 IF(IFSET)891,895,895
  891 DO 892 I = 1,NX
  892 A(2*I) = -A(2*I)
  895 RETURN
C
C     THE FOLLOWING PROGRAM COMPUTES THE SIN AND INV TABLES.
C
  900 MT=MAX0(M(1),M(2),M(3)) -2
      MT = MAX0(2,MT)
  904 IF (MT-18) 906,906,13
  906 IFERR=0
      NT=2**MT
      NTV2=NT/2
C
C     SET UP SIN TABLE
C     THETA=PIE/2**(L+1) FOR L=1
  910 THETA=.7853981634
C
C     JSTEP=2**(MT-L+1) FOR L=1
      JSTEP=NT
C
C     JDIF=2**(MT-L) FOR L=1
      JDIF=NTV2
      S(JDIF)=SIN(THETA)
      DO 950 L=2,MT
      THETA=THETA/2.
      JSTEP2=JSTEP
      JSTEP=JDIF
      JDIF=JSTEP/2
      S(JDIF)=SIN(THETA)
      JC1=NT-JDIF
      S(JC1)=COS(THETA)
      JLAST=NT-JSTEP2
      IF(JLAST - JSTEP) 950,920,920
  920 DO 940 J=JSTEP,JLAST,JSTEP
      JC=NT-J
      JD=J+JDIF
  940 S(JD)=S(J)*S(JC1)+S(JDIF)*S(JC)
  950 CONTINUE
C
C     SET UP INV(J) TABLE
C
  960 MTLEXP=NTV2
C
C     MTLEXP=2**(MT-L). FOR L=1
      LM1EXP=1
C
C     LM1EXP=2**(L-1). FOR L=1
      INV(1)=0
      DO 980 L=1,MT
      INV(LM1EXP+1) = MTLEXP
      DO 970 J=2,LM1EXP
      JJ=J+LM1EXP
  970 INV(JJ)=INV(J)+MTLEXP
      MTLEXP=MTLEXP/2
  980 LM1EXP=LM1EXP*2
  982 IF(IFSET)12,895,12
      END