1***************************** -*- Mode: Fortran -*- ********************
2*
3* cfft.f -- FFT routines.
4*
5* Copyright (C) 1994-95  K. Scott Hunziker.
6* Copyright (C) 1990-94  The Boeing Company.
7*
8* See the file COPYING for license, warranty, and permission details.
9*
10************************************************************************
11*
12* $Id: cfft.f,v 1.2 1996/10/19 19:55:37 ksh Exp $
13*
14      SUBROUTINE CFFTI(N,WSAVE)
15CSE
16      IMPLICIT REAL*8(A-H,O-Z)
17C***BEGIN PROLOGUE  CFFTI
18C***DATE WRITTEN   790601   (YYMMDD)
19C***REVISION DATE  830401   (YYMMDD)
20C***CATEGORY NO.  J1A2
21C***KEYWORDS  FOURIER TRANSFORM
22C***AUTHOR  SWARZTRAUBER, P. N., (NCAR)
23CPS
24C***PURPOSE  Initialize for CFFTF and CFFTB.
25C***DESCRIPTION
26C
27C  Subroutine CFFTI initializes the array WSAVE which is used in
28C  both CFFTF and CFFTB.  The prime factorization of N together with
29C  a tabulation of the trigonometric functions are computed and
30C  stored in WSAVE.
31CPE
32CAS
33C  Input Parameter
34C
35C  N       the length of the sequence to be transformed
36C
37C  Output Parameter
38C
39C  WSAVE   a work array which must be dimensioned at least 4*N+15.
40C          The same work array can be used for both CFFTF and CFFTB
41C          as long as N remains unchanged.  Different WSAVE arrays
42C          are required for different values of N.  The contents of
43C          WSAVE must not be changed between calls of CFFTF or CFFTB.
44CAE
45C***REFERENCES  (NONE)
46C***ROUTINES CALLED  CFFTI1
47C***END PROLOGUE  CFFTI
48      DIMENSION       WSAVE(1)
49C***FIRST EXECUTABLE STATEMENT  CFFTI
50      IF (N .EQ. 1) RETURN
51      IW1 = N+N+1
52      IW2 = IW1+N+N
53      CALL CFFTI1 (N,WSAVE(IW1),WSAVE(IW2))
54      RETURN
55      END
56*
57      SUBROUTINE CFFTI1(N,WA,IFAC)
58CSE
59      IMPLICIT REAL*8(A-H,O-Z)
60C***BEGIN PROLOGUE  CFFTI1
61C***REFER TO  CFFTI
62C***ROUTINES CALLED  (NONE)
63C***END PROLOGUE  CFFTI1
64      DIMENSION       WA(1)      ,IFAC(3)    ,NTRYH(4)
65      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/3,4,2,5/
66C***FIRST EXECUTABLE STATEMENT  CFFTI1
67      NL = N
68      NF = 0
69      J = 0
70  101 J = J+1
71      IF (J-4) 102,102,103
72  102 NTRY = NTRYH(J)
73      GO TO 104
74  103 NTRY = NTRY+2
75  104 NQ = NL/NTRY
76      NR = NL-NTRY*NQ
77      IF (NR) 101,105,101
78  105 NF = NF+1
79      IFAC(NF+2) = NTRY
80      NL = NQ
81      IF (NTRY .NE. 2) GO TO 107
82      IF (NF .EQ. 1) GO TO 107
83      DO 106 I=2,NF
84         IB = NF-I+2
85         IFAC(IB+2) = IFAC(IB+1)
86  106 CONTINUE
87      IFAC(3) = 2
88  107 IF (NL .NE. 1) GO TO 104
89      IFAC(1) = N
90      IFAC(2) = NF
91      TPI = 6.28318530717959
92      ARGH = TPI/DBLE(N)
93      I = 2
94      L1 = 1
95      DO 110 K1=1,NF
96         IP = IFAC(K1+2)
97         LD = 0
98         L2 = L1*IP
99         IDO = N/L2
100         IDOT = IDO+IDO+2
101         IPM = IP-1
102         DO 109 J=1,IPM
103            I1 = I
104            WA(I-1) = 1.
105            WA(I) = 0.
106            LD = LD+L1
107            FI = 0.
108            ARGLD = DBLE(LD)*ARGH
109            DO 108 II=4,IDOT,2
110               I = I+2
111               FI = FI+1.
112               ARG = FI*ARGLD
113               WA(I-1) = COS(ARG)
114               WA(I) = SIN(ARG)
115  108       CONTINUE
116            IF (IP .LE. 5) GO TO 109
117            WA(I1-1) = WA(I-1)
118            WA(I1) = WA(I)
119  109    CONTINUE
120         L1 = L2
121  110 CONTINUE
122      RETURN
123      END
124*
125      SUBROUTINE CFFTF(N,C,WSAVE)
126CSE
127      IMPLICIT REAL*8(A-H,O-Z)
128C***BEGIN PROLOGUE  CFFTF
129C***DATE WRITTEN   790601   (YYMMDD)
130C***REVISION DATE  800626   (YYMMDD)
131C***CATEGORY NO.  J1A2
132C***KEYWORDS  FOURIER TRANSFORM
133C***AUTHOR  SWARZTRAUBER, P. N., (NCAR)
134CPS
135C***PURPOSE  Forward transform of a complex, periodic sequence.
136C***DESCRIPTION
137C
138C  Subroutine CFFTF computes the forward complex discrete Fourier
139C  transform (the Fourier analysis).  Equivalently, CFFTF computes
140C  the Fourier coefficients of a complex periodic sequence.
141C  The transform is defined below at output parameter C.
142C
143C  The transform is not normalized.  To obtain a normalized transform
144C  the output must be divided by N.  Otherwise a call of CFFTF
145C  followed by a call of CFFTB will multiply the sequence by N.
146C
147C  The array WSAVE which is used by subroutine CFFTF must be
148C  initialized by calling subroutine CFFTI(N,WSAVE).
149CPE
150CAS
151C  Input Parameters
152C
153C
154C  N      the length of the complex sequence C.  The method is
155C         more efficient when N is the product of small primes.
156C
157C  C      a complex array of length N which contains the sequence
158C
159C  WSAVE   a real work array which must be dimensioned at least 4*N+15
160C          in the program that calls CFFTF.  The WSAVE array must be
161C          initialized by calling subroutine CFFTI(N,WSAVE), and a
162C          different WSAVE array must be used for each different
163C          value of N.  This initialization does not have to be
164C          repeated so long as N remains unchanged.  Thus subsequent
165C          transforms can be obtained faster than the first.
166C          The same WSAVE array can be used by CFFTF and CFFTB.
167C
168C  Output Parameters
169C
170C  C      for J=1,...,N
171C
172C             C(J)=the sum from K=1,...,N of
173C
174C                   C(K)*EXP(-I*J*K*2*PI/N)
175C
176C                         where I=SQRT(-1)
177C
178C  WSAVE   contains initialization calculations which must not be
179C          destroyed between calls of subroutine CFFTF or CFFTB
180CAE
181C***REFERENCES  (NONE)
182C***ROUTINES CALLED  CFFTF1
183C***END PROLOGUE  CFFTF
184      DIMENSION       C(1)       ,WSAVE(1)
185C***FIRST EXECUTABLE STATEMENT  CFFTF
186      IF (N .EQ. 1) RETURN
187      IW1 = N+N+1
188      IW2 = IW1+N+N
189      CALL CFFTF1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
190      RETURN
191      END
192*
193      SUBROUTINE CFFTF1(N,C,CH,WA,IFAC)
194CSE
195      IMPLICIT REAL*8(A-H,O-Z)
196C***BEGIN PROLOGUE  CFFTF1
197C***REFER TO  CFFTF
198C***ROUTINES CALLED  PASSF,PASSF2,PASSF3,PASSF4,PASSF5
199C***END PROLOGUE  CFFTF1
200      DIMENSION       CH(1)      ,C(1)       ,WA(1)      ,IFAC(2)
201C***FIRST EXECUTABLE STATEMENT  CFFTF1
202      NF = IFAC(2)
203      NA = 0
204      L1 = 1
205      IW = 1
206      DO 116 K1=1,NF
207         IP = IFAC(K1+2)
208         L2 = IP*L1
209         IDO = N/L2
210         IDOT = IDO+IDO
211         IDL1 = IDOT*L1
212         IF (IP .NE. 4) GO TO 103
213         IX2 = IW+IDOT
214         IX3 = IX2+IDOT
215         IF (NA .NE. 0) GO TO 101
216         CALL PASSF4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
217         GO TO 102
218  101    CALL PASSF4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
219  102    NA = 1-NA
220         GO TO 115
221  103    IF (IP .NE. 2) GO TO 106
222         IF (NA .NE. 0) GO TO 104
223         CALL PASSF2 (IDOT,L1,C,CH,WA(IW))
224         GO TO 105
225  104    CALL PASSF2 (IDOT,L1,CH,C,WA(IW))
226  105    NA = 1-NA
227         GO TO 115
228  106    IF (IP .NE. 3) GO TO 109
229         IX2 = IW+IDOT
230         IF (NA .NE. 0) GO TO 107
231         CALL PASSF3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
232         GO TO 108
233  107    CALL PASSF3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
234  108    NA = 1-NA
235         GO TO 115
236  109    IF (IP .NE. 5) GO TO 112
237         IX2 = IW+IDOT
238         IX3 = IX2+IDOT
239         IX4 = IX3+IDOT
240         IF (NA .NE. 0) GO TO 110
241         CALL PASSF5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
242         GO TO 111
243  110    CALL PASSF5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
244  111    NA = 1-NA
245         GO TO 115
246  112    IF (NA .NE. 0) GO TO 113
247         CALL PASSF (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
248         GO TO 114
249  113    CALL PASSF (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
250  114    IF (NAC .NE. 0) NA = 1-NA
251  115    L1 = L2
252         IW = IW+(IP-1)*IDOT
253  116 CONTINUE
254      IF (NA .EQ. 0) RETURN
255      N2 = N+N
256      DO 117 I=1,N2
257         C(I) = CH(I)
258  117 CONTINUE
259      RETURN
260      END
261*
262      SUBROUTINE PASSF(NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
263CSE
264      IMPLICIT REAL*8(A-H,O-Z)
265C***BEGIN PROLOGUE  PASSF
266C***REFER TO  CFFTF
267C***ROUTINES CALLED  (NONE)
268C***END PROLOGUE  PASSF
269      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
270     1                C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP),
271     2                CH2(IDL1,IP)
272C***FIRST EXECUTABLE STATEMENT  PASSF
273      IDOT = IDO/2
274      NT = IP*IDL1
275      IPP2 = IP+2
276      IPPH = (IP+1)/2
277      IDP = IP*IDO
278C
279      IF (IDO .LT. L1) GO TO 106
280      DO 103 J=2,IPPH
281         JC = IPP2-J
282         DO 102 K=1,L1
283CDIR$ IVDEP
284            DO 101 I=1,IDO
285               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
286               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
287  101       CONTINUE
288  102    CONTINUE
289  103 CONTINUE
290      DO 105 K=1,L1
291CDIR$ IVDEP
292         DO 104 I=1,IDO
293            CH(I,K,1) = CC(I,1,K)
294  104    CONTINUE
295  105 CONTINUE
296      GO TO 112
297  106 DO 109 J=2,IPPH
298         JC = IPP2-J
299         DO 108 I=1,IDO
300CDIR$ IVDEP
301            DO 107 K=1,L1
302               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
303               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
304  107       CONTINUE
305  108    CONTINUE
306  109 CONTINUE
307      DO 111 I=1,IDO
308CDIR$ IVDEP
309         DO 110 K=1,L1
310            CH(I,K,1) = CC(I,1,K)
311  110    CONTINUE
312  111 CONTINUE
313  112 IDL = 2-IDO
314      INC = 0
315      DO 116 L=2,IPPH
316         LC = IPP2-L
317         IDL = IDL+IDO
318CDIR$ IVDEP
319         DO 113 IK=1,IDL1
320            C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
321            C2(IK,LC) = -WA(IDL)*CH2(IK,IP)
322  113    CONTINUE
323         IDLJ = IDL
324         INC = INC+IDO
325         DO 115 J=3,IPPH
326            JC = IPP2-J
327            IDLJ = IDLJ+INC
328            IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP
329            WAR = WA(IDLJ-1)
330            WAI = WA(IDLJ)
331CDIR$ IVDEP
332            DO 114 IK=1,IDL1
333               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
334               C2(IK,LC) = C2(IK,LC)-WAI*CH2(IK,JC)
335  114       CONTINUE
336  115    CONTINUE
337  116 CONTINUE
338      DO 118 J=2,IPPH
339CDIR$ IVDEP
340         DO 117 IK=1,IDL1
341            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
342  117    CONTINUE
343  118 CONTINUE
344      DO 120 J=2,IPPH
345         JC = IPP2-J
346CDIR$ IVDEP
347         DO 119 IK=2,IDL1,2
348            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
349            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
350            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
351            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
352  119    CONTINUE
353  120 CONTINUE
354      NAC = 1
355      IF (IDO .EQ. 2) RETURN
356      NAC = 0
357CDIR$ IVDEP
358      DO 121 IK=1,IDL1
359         C2(IK,1) = CH2(IK,1)
360  121 CONTINUE
361      DO 123 J=2,IP
362CDIR$ IVDEP
363         DO 122 K=1,L1
364            C1(1,K,J) = CH(1,K,J)
365            C1(2,K,J) = CH(2,K,J)
366  122    CONTINUE
367  123 CONTINUE
368      IF (IDOT .GT. L1) GO TO 127
369      IDIJ = 0
370      DO 126 J=2,IP
371         IDIJ = IDIJ+2
372         DO 125 I=4,IDO,2
373            IDIJ = IDIJ+2
374CDIR$ IVDEP
375            DO 124 K=1,L1
376               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
377               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
378  124       CONTINUE
379  125    CONTINUE
380  126 CONTINUE
381      RETURN
382  127 IDJ = 2-IDO
383      DO 130 J=2,IP
384         IDJ = IDJ+IDO
385         DO 129 K=1,L1
386            IDIJ = IDJ
387CDIR$ IVDEP
388            DO 128 I=4,IDO,2
389               IDIJ = IDIJ+2
390               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)+WA(IDIJ)*CH(I,K,J)
391               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)-WA(IDIJ)*CH(I-1,K,J)
392  128       CONTINUE
393  129    CONTINUE
394  130 CONTINUE
395      RETURN
396      END
397*
398      SUBROUTINE PASSF2(IDO,L1,CC,CH,WA1)
399CSE
400      IMPLICIT REAL*8(A-H,O-Z)
401C***BEGIN PROLOGUE  PASSF2
402C***REFER TO  CFFTF
403C***ROUTINES CALLED  (NONE)
404C***END PROLOGUE  PASSF2
405      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)           ,
406     1                WA1(1)
407C***FIRST EXECUTABLE STATEMENT  PASSF2
408      IF (IDO .GT. 2) GO TO 102
409      DO 101 K=1,L1
410         CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
411         CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
412         CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
413         CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
414  101 CONTINUE
415      RETURN
416  102 IF(IDO/2.LT.L1) GO TO 105
417      DO 104 K=1,L1
418CDIR$ IVDEP
419         DO 103 I=2,IDO,2
420            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
421            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
422            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
423            TI2 = CC(I,1,K)-CC(I,2,K)
424            CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2
425            CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2
426  103    CONTINUE
427  104 CONTINUE
428      RETURN
429  105 DO 107 I=2,IDO,2
430CDIR$ IVDEP
431      DO 106 K=1,L1
432            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
433            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
434            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
435            TI2 = CC(I,1,K)-CC(I,2,K)
436            CH(I,K,2) = WA1(I-1)*TI2-WA1(I)*TR2
437            CH(I-1,K,2) = WA1(I-1)*TR2+WA1(I)*TI2
438  106    CONTINUE
439  107 CONTINUE
440      RETURN
441      END
442*
443      SUBROUTINE PASSF3(IDO,L1,CC,CH,WA1,WA2)
444CSE
445      IMPLICIT REAL*8(A-H,O-Z)
446C***BEGIN PROLOGUE  PASSF3
447C***REFER TO  CFFTF
448C***ROUTINES CALLED  (NONE)
449C***END PROLOGUE  PASSF3
450      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)           ,
451     1                WA1(1)     ,WA2(1)
452      DATA TAUR,TAUI /-.5,-.866025403784439/
453C***FIRST EXECUTABLE STATEMENT  PASSF3
454      IF (IDO .NE. 2) GO TO 102
455      DO 101 K=1,L1
456         TR2 = CC(1,2,K)+CC(1,3,K)
457         CR2 = CC(1,1,K)+TAUR*TR2
458         CH(1,K,1) = CC(1,1,K)+TR2
459         TI2 = CC(2,2,K)+CC(2,3,K)
460         CI2 = CC(2,1,K)+TAUR*TI2
461         CH(2,K,1) = CC(2,1,K)+TI2
462         CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
463         CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
464         CH(1,K,2) = CR2-CI3
465         CH(1,K,3) = CR2+CI3
466         CH(2,K,2) = CI2+CR3
467         CH(2,K,3) = CI2-CR3
468  101 CONTINUE
469      RETURN
470  102 IF(IDO/2.LT.L1) GO TO 105
471      DO 104 K=1,L1
472CDIR$ IVDEP
473         DO 103 I=2,IDO,2
474            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
475            CR2 = CC(I-1,1,K)+TAUR*TR2
476            CH(I-1,K,1) = CC(I-1,1,K)+TR2
477            TI2 = CC(I,2,K)+CC(I,3,K)
478            CI2 = CC(I,1,K)+TAUR*TI2
479            CH(I,K,1) = CC(I,1,K)+TI2
480            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
481            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
482            DR2 = CR2-CI3
483            DR3 = CR2+CI3
484            DI2 = CI2+CR3
485            DI3 = CI2-CR3
486            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
487            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
488            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
489            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
490  103    CONTINUE
491  104 CONTINUE
492      RETURN
493  105 DO 107 I=2,IDO,2
494CDIR$ IVDEP
495         DO 106 K=1,L1
496            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
497            CR2 = CC(I-1,1,K)+TAUR*TR2
498            CH(I-1,K,1) = CC(I-1,1,K)+TR2
499            TI2 = CC(I,2,K)+CC(I,3,K)
500            CI2 = CC(I,1,K)+TAUR*TI2
501            CH(I,K,1) = CC(I,1,K)+TI2
502            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
503            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
504            DR2 = CR2-CI3
505            DR3 = CR2+CI3
506            DI2 = CI2+CR3
507            DI3 = CI2-CR3
508            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
509            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
510            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
511            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
512  106    CONTINUE
513  107 CONTINUE
514      RETURN
515      END
516*
517      SUBROUTINE PASSF4(IDO,L1,CC,CH,WA1,WA2,WA3)
518CSE
519      IMPLICIT REAL*8(A-H,O-Z)
520C***BEGIN PROLOGUE  PASSF4
521C***REFER TO  CFFTF
522C***ROUTINES CALLED  (NONE)
523C***END PROLOGUE  PASSF4
524      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)           ,
525     1                WA1(1)     ,WA2(1)     ,WA3(1)
526C***FIRST EXECUTABLE STATEMENT  PASSF4
527      IF (IDO .NE. 2) GO TO 102
528      DO 101 K=1,L1
529         TI1 = CC(2,1,K)-CC(2,3,K)
530         TI2 = CC(2,1,K)+CC(2,3,K)
531         TR4 = CC(2,2,K)-CC(2,4,K)
532         TI3 = CC(2,2,K)+CC(2,4,K)
533         TR1 = CC(1,1,K)-CC(1,3,K)
534         TR2 = CC(1,1,K)+CC(1,3,K)
535         TI4 = CC(1,4,K)-CC(1,2,K)
536         TR3 = CC(1,2,K)+CC(1,4,K)
537         CH(1,K,1) = TR2+TR3
538         CH(1,K,3) = TR2-TR3
539         CH(2,K,1) = TI2+TI3
540         CH(2,K,3) = TI2-TI3
541         CH(1,K,2) = TR1+TR4
542         CH(1,K,4) = TR1-TR4
543         CH(2,K,2) = TI1+TI4
544         CH(2,K,4) = TI1-TI4
545  101 CONTINUE
546      RETURN
547  102 IF(IDO/2.LT.L1) GO TO 105
548      DO 104 K=1,L1
549CDIR$ IVDEP
550         DO 103 I=2,IDO,2
551            TI1 = CC(I,1,K)-CC(I,3,K)
552            TI2 = CC(I,1,K)+CC(I,3,K)
553            TI3 = CC(I,2,K)+CC(I,4,K)
554            TR4 = CC(I,2,K)-CC(I,4,K)
555            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
556            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
557            TI4 = CC(I-1,4,K)-CC(I-1,2,K)
558            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
559            CH(I-1,K,1) = TR2+TR3
560            CR3 = TR2-TR3
561            CH(I,K,1) = TI2+TI3
562            CI3 = TI2-TI3
563            CR2 = TR1+TR4
564            CR4 = TR1-TR4
565            CI2 = TI1+TI4
566            CI4 = TI1-TI4
567            CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2
568            CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2
569            CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3
570            CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3
571            CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4
572            CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4
573  103    CONTINUE
574  104 CONTINUE
575      RETURN
576  105 DO 107 I=2,IDO,2
577CDIR$ IVDEP
578         DO 106 K=1,L1
579            TI1 = CC(I,1,K)-CC(I,3,K)
580            TI2 = CC(I,1,K)+CC(I,3,K)
581            TI3 = CC(I,2,K)+CC(I,4,K)
582            TR4 = CC(I,2,K)-CC(I,4,K)
583            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
584            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
585            TI4 = CC(I-1,4,K)-CC(I-1,2,K)
586            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
587            CH(I-1,K,1) = TR2+TR3
588            CR3 = TR2-TR3
589            CH(I,K,1) = TI2+TI3
590            CI3 = TI2-TI3
591            CR2 = TR1+TR4
592            CR4 = TR1-TR4
593            CI2 = TI1+TI4
594            CI4 = TI1-TI4
595            CH(I-1,K,2) = WA1(I-1)*CR2+WA1(I)*CI2
596            CH(I,K,2) = WA1(I-1)*CI2-WA1(I)*CR2
597            CH(I-1,K,3) = WA2(I-1)*CR3+WA2(I)*CI3
598            CH(I,K,3) = WA2(I-1)*CI3-WA2(I)*CR3
599            CH(I-1,K,4) = WA3(I-1)*CR4+WA3(I)*CI4
600            CH(I,K,4) = WA3(I-1)*CI4-WA3(I)*CR4
601  106    CONTINUE
602  107 CONTINUE
603      RETURN
604      END
605*
606      SUBROUTINE PASSF5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
607CSE
608      IMPLICIT REAL*8(A-H,O-Z)
609C***BEGIN PROLOGUE  PASSF5
610C***REFER TO  CFFTF
611C***ROUTINES CALLED  (NONE)
612C***END PROLOGUE  PASSF5
613      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)           ,
614     1                WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
615      DATA TR11,TI11,TR12,TI12 /.309016994374947,-.951056516295154,
616     1-.809016994374947,-.587785252292473/
617C***FIRST EXECUTABLE STATEMENT  PASSF5
618      IF (IDO .NE. 2) GO TO 102
619      DO 101 K=1,L1
620         TI5 = CC(2,2,K)-CC(2,5,K)
621         TI2 = CC(2,2,K)+CC(2,5,K)
622         TI4 = CC(2,3,K)-CC(2,4,K)
623         TI3 = CC(2,3,K)+CC(2,4,K)
624         TR5 = CC(1,2,K)-CC(1,5,K)
625         TR2 = CC(1,2,K)+CC(1,5,K)
626         TR4 = CC(1,3,K)-CC(1,4,K)
627         TR3 = CC(1,3,K)+CC(1,4,K)
628         CH(1,K,1) = CC(1,1,K)+TR2+TR3
629         CH(2,K,1) = CC(2,1,K)+TI2+TI3
630         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
631         CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
632         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
633         CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
634         CR5 = TI11*TR5+TI12*TR4
635         CI5 = TI11*TI5+TI12*TI4
636         CR4 = TI12*TR5-TI11*TR4
637         CI4 = TI12*TI5-TI11*TI4
638         CH(1,K,2) = CR2-CI5
639         CH(1,K,5) = CR2+CI5
640         CH(2,K,2) = CI2+CR5
641         CH(2,K,3) = CI3+CR4
642         CH(1,K,3) = CR3-CI4
643         CH(1,K,4) = CR3+CI4
644         CH(2,K,4) = CI3-CR4
645         CH(2,K,5) = CI2-CR5
646  101 CONTINUE
647      RETURN
648  102 IF(IDO/2.LT.L1) GO TO 105
649      DO 104 K=1,L1
650CDIR$ IVDEP
651         DO 103 I=2,IDO,2
652            TI5 = CC(I,2,K)-CC(I,5,K)
653            TI2 = CC(I,2,K)+CC(I,5,K)
654            TI4 = CC(I,3,K)-CC(I,4,K)
655            TI3 = CC(I,3,K)+CC(I,4,K)
656            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
657            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
658            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
659            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
660            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
661            CH(I,K,1) = CC(I,1,K)+TI2+TI3
662            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
663            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
664            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
665            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
666            CR5 = TI11*TR5+TI12*TR4
667            CI5 = TI11*TI5+TI12*TI4
668            CR4 = TI12*TR5-TI11*TR4
669            CI4 = TI12*TI5-TI11*TI4
670            DR3 = CR3-CI4
671            DR4 = CR3+CI4
672            DI3 = CI3+CR4
673            DI4 = CI3-CR4
674            DR5 = CR2+CI5
675            DR2 = CR2-CI5
676            DI5 = CI2-CR5
677            DI2 = CI2+CR5
678            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
679            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
680            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
681            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
682            CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4
683            CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4
684            CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5
685            CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5
686  103    CONTINUE
687  104 CONTINUE
688      RETURN
689  105 DO 107 I=2,IDO,2
690CDIR$ IVDEP
691         DO 106 K=1,L1
692            TI5 = CC(I,2,K)-CC(I,5,K)
693            TI2 = CC(I,2,K)+CC(I,5,K)
694            TI4 = CC(I,3,K)-CC(I,4,K)
695            TI3 = CC(I,3,K)+CC(I,4,K)
696            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
697            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
698            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
699            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
700            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
701            CH(I,K,1) = CC(I,1,K)+TI2+TI3
702            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
703            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
704            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
705            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
706            CR5 = TI11*TR5+TI12*TR4
707            CI5 = TI11*TI5+TI12*TI4
708            CR4 = TI12*TR5-TI11*TR4
709            CI4 = TI12*TI5-TI11*TI4
710            DR3 = CR3-CI4
711            DR4 = CR3+CI4
712            DI3 = CI3+CR4
713            DI4 = CI3-CR4
714            DR5 = CR2+CI5
715            DR2 = CR2-CI5
716            DI5 = CI2-CR5
717            DI2 = CI2+CR5
718            CH(I-1,K,2) = WA1(I-1)*DR2+WA1(I)*DI2
719            CH(I,K,2) = WA1(I-1)*DI2-WA1(I)*DR2
720            CH(I-1,K,3) = WA2(I-1)*DR3+WA2(I)*DI3
721            CH(I,K,3) = WA2(I-1)*DI3-WA2(I)*DR3
722            CH(I-1,K,4) = WA3(I-1)*DR4+WA3(I)*DI4
723            CH(I,K,4) = WA3(I-1)*DI4-WA3(I)*DR4
724            CH(I-1,K,5) = WA4(I-1)*DR5+WA4(I)*DI5
725            CH(I,K,5) = WA4(I-1)*DI5-WA4(I)*DR5
726  106    CONTINUE
727  107 CONTINUE
728      RETURN
729      END
730*
731      SUBROUTINE CFFTB(N,C,WSAVE)
732CSE
733      IMPLICIT REAL*8(A-H,O-Z)
734C***BEGIN PROLOGUE  CFFTB
735C***DATE WRITTEN   790601   (YYMMDD)
736C***REVISION DATE  830401   (YYMMDD)
737C***CATEGORY NO.  J1A2
738C***KEYWORDS  FOURIER TRANSFORM
739C***AUTHOR  SWARZTRAUBER, P. N., (NCAR)
740CPS
741C***PURPOSE  Unnormalized inverse of CFFTF.
742C***DESCRIPTION
743C
744C  Subroutine CFFTB computes the backward complex discrete Fourier
745C  transform (the Fourier synthesis).  Equivalently, CFFTB computes
746C  a complex periodic sequence from its Fourier coefficients.
747C  The transform is defined below at output parameter C.
748C
749C  A call of CFFTF followed by a call of CFFTB will multiply the
750C  sequence by N.
751C
752C  The array WSAVE which is used by subroutine CFFTB must be
753C  initialized by calling subroutine CFFTI(N,WSAVE).
754CPE
755CAS
756C  Input Parameters
757C
758C
759C  N      the length of the complex sequence C.  The method is
760C         more efficient when N is the product of small primes.
761C
762C  C      a complex array of length N which contains the sequence
763C
764C  WSAVE   a real work array which must be dimensioned at least 4*N+15
765C          in the program that calls CFFTB.  The WSAVE array must be
766C          initialized by calling subroutine CFFTI(N,WSAVE), and a
767C          different WSAVE array must be used for each different
768C          value of N.  This initialization does not have to be
769C          repeated so long as N remains unchanged.  Thus subsequent
770C          transforms can be obtained faster than the first.
771C          The same WSAVE array can be used by CFFTF and CFFTB.
772C
773C  Output Parameters
774C
775C  C      For J=1,...,N
776C
777C             C(J)=the sum from K=1,...,N of
778C
779C                   C(K)*EXP(I*J*K*2*PI/N)
780C
781C                         where I=SQRT(-1)
782C
783C  WSAVE   contains initialization calculations which must not be
784C          destroyed between calls of subroutine CFFTF or CFFTB
785CAE
786C***REFERENCES  (NONE)
787C***ROUTINES CALLED  CFFTB1
788C***END PROLOGUE  CFFTB
789      DIMENSION       C(1)       ,WSAVE(1)
790C***FIRST EXECUTABLE STATEMENT  CFFTB
791      IF (N .EQ. 1) RETURN
792      IW1 = N+N+1
793      IW2 = IW1+N+N
794      CALL CFFTB1 (N,C,WSAVE,WSAVE(IW1),WSAVE(IW2))
795      RETURN
796      END
797*
798      SUBROUTINE CFFTB1(N,C,CH,WA,IFAC)
799CSE
800      IMPLICIT REAL*8(A-H,O-Z)
801C***BEGIN PROLOGUE  CFFTB1
802C***REFER TO  CFFTB
803C***ROUTINES CALLED  PASSB,PASSB2,PASSB3,PASSB4,PASSB5
804C***END PROLOGUE  CFFTB1
805      DIMENSION       CH(1)      ,C(1)       ,WA(1)      ,IFAC(2)
806C***FIRST EXECUTABLE STATEMENT  CFFTB1
807      NF = IFAC(2)
808      NA = 0
809      L1 = 1
810      IW = 1
811      DO 116 K1=1,NF
812         IP = IFAC(K1+2)
813         L2 = IP*L1
814         IDO = N/L2
815         IDOT = IDO+IDO
816         IDL1 = IDOT*L1
817         IF (IP .NE. 4) GO TO 103
818         IX2 = IW+IDOT
819         IX3 = IX2+IDOT
820         IF (NA .NE. 0) GO TO 101
821         CALL PASSB4 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
822         GO TO 102
823  101    CALL PASSB4 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
824  102    NA = 1-NA
825         GO TO 115
826  103    IF (IP .NE. 2) GO TO 106
827         IF (NA .NE. 0) GO TO 104
828         CALL PASSB2 (IDOT,L1,C,CH,WA(IW))
829         GO TO 105
830  104    CALL PASSB2 (IDOT,L1,CH,C,WA(IW))
831  105    NA = 1-NA
832         GO TO 115
833  106    IF (IP .NE. 3) GO TO 109
834         IX2 = IW+IDOT
835         IF (NA .NE. 0) GO TO 107
836         CALL PASSB3 (IDOT,L1,C,CH,WA(IW),WA(IX2))
837         GO TO 108
838  107    CALL PASSB3 (IDOT,L1,CH,C,WA(IW),WA(IX2))
839  108    NA = 1-NA
840         GO TO 115
841  109    IF (IP .NE. 5) GO TO 112
842         IX2 = IW+IDOT
843         IX3 = IX2+IDOT
844         IX4 = IX3+IDOT
845         IF (NA .NE. 0) GO TO 110
846         CALL PASSB5 (IDOT,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
847         GO TO 111
848  110    CALL PASSB5 (IDOT,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
849  111    NA = 1-NA
850         GO TO 115
851  112    IF (NA .NE. 0) GO TO 113
852         CALL PASSB (NAC,IDOT,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
853         GO TO 114
854  113    CALL PASSB (NAC,IDOT,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
855  114    IF (NAC .NE. 0) NA = 1-NA
856  115    L1 = L2
857         IW = IW+(IP-1)*IDOT
858  116 CONTINUE
859      IF (NA .EQ. 0) RETURN
860      N2 = N+N
861      DO 117 I=1,N2
862         C(I) = CH(I)
863  117 CONTINUE
864      RETURN
865      END
866*
867      SUBROUTINE PASSB(NAC,IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
868CSE
869      IMPLICIT REAL*8(A-H,O-Z)
870C***BEGIN PROLOGUE  PASSB
871C***REFER TO  CFFTB
872C***ROUTINES CALLED  (NONE)
873C***END PROLOGUE  PASSB
874      DIMENSION       CH(IDO,L1,IP)          ,CC(IDO,IP,L1)          ,
875     1                C1(IDO,L1,IP)          ,WA(1)      ,C2(IDL1,IP),
876     2                CH2(IDL1,IP)
877C***FIRST EXECUTABLE STATEMENT  PASSB
878      IDOT = IDO/2
879      NT = IP*IDL1
880      IPP2 = IP+2
881      IPPH = (IP+1)/2
882      IDP = IP*IDO
883C
884      IF (IDO .LT. L1) GO TO 106
885      DO 103 J=2,IPPH
886         JC = IPP2-J
887         DO 102 K=1,L1
888CDIR$ IVDEP
889            DO 101 I=1,IDO
890               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
891               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
892  101       CONTINUE
893  102    CONTINUE
894  103 CONTINUE
895      DO 105 K=1,L1
896CDIR$ IVDEP
897         DO 104 I=1,IDO
898            CH(I,K,1) = CC(I,1,K)
899  104    CONTINUE
900  105 CONTINUE
901      GO TO 112
902  106 DO 109 J=2,IPPH
903         JC = IPP2-J
904         DO 108 I=1,IDO
905CDIR$ IVDEP
906            DO 107 K=1,L1
907               CH(I,K,J) = CC(I,J,K)+CC(I,JC,K)
908               CH(I,K,JC) = CC(I,J,K)-CC(I,JC,K)
909  107       CONTINUE
910  108    CONTINUE
911  109 CONTINUE
912      DO 111 I=1,IDO
913CDIR$ IVDEP
914         DO 110 K=1,L1
915            CH(I,K,1) = CC(I,1,K)
916  110    CONTINUE
917  111 CONTINUE
918  112 IDL = 2-IDO
919      INC = 0
920      DO 116 L=2,IPPH
921         LC = IPP2-L
922         IDL = IDL+IDO
923CDIR$ IVDEP
924         DO 113 IK=1,IDL1
925            C2(IK,L) = CH2(IK,1)+WA(IDL-1)*CH2(IK,2)
926            C2(IK,LC) = WA(IDL)*CH2(IK,IP)
927  113    CONTINUE
928         IDLJ = IDL
929         INC = INC+IDO
930         DO 115 J=3,IPPH
931            JC = IPP2-J
932            IDLJ = IDLJ+INC
933            IF (IDLJ .GT. IDP) IDLJ = IDLJ-IDP
934            WAR = WA(IDLJ-1)
935            WAI = WA(IDLJ)
936CDIR$ IVDEP
937            DO 114 IK=1,IDL1
938               C2(IK,L) = C2(IK,L)+WAR*CH2(IK,J)
939               C2(IK,LC) = C2(IK,LC)+WAI*CH2(IK,JC)
940  114       CONTINUE
941  115    CONTINUE
942  116 CONTINUE
943      DO 118 J=2,IPPH
944CDIR$ IVDEP
945         DO 117 IK=1,IDL1
946            CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
947  117    CONTINUE
948  118 CONTINUE
949      DO 120 J=2,IPPH
950         JC = IPP2-J
951CDIR$ IVDEP
952         DO 119 IK=2,IDL1,2
953            CH2(IK-1,J) = C2(IK-1,J)-C2(IK,JC)
954            CH2(IK-1,JC) = C2(IK-1,J)+C2(IK,JC)
955            CH2(IK,J) = C2(IK,J)+C2(IK-1,JC)
956            CH2(IK,JC) = C2(IK,J)-C2(IK-1,JC)
957  119    CONTINUE
958  120 CONTINUE
959      NAC = 1
960      IF (IDO .EQ. 2) RETURN
961      NAC = 0
962      DO 121 IK=1,IDL1
963         C2(IK,1) = CH2(IK,1)
964  121 CONTINUE
965      DO 123 J=2,IP
966CDIR$ IVDEP
967         DO 122 K=1,L1
968            C1(1,K,J) = CH(1,K,J)
969            C1(2,K,J) = CH(2,K,J)
970  122    CONTINUE
971  123 CONTINUE
972      IF (IDOT .GT. L1) GO TO 127
973      IDIJ = 0
974      DO 126 J=2,IP
975         IDIJ = IDIJ+2
976         DO 125 I=4,IDO,2
977            IDIJ = IDIJ+2
978CDIR$ IVDEP
979            DO 124 K=1,L1
980               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
981               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
982  124       CONTINUE
983  125    CONTINUE
984  126 CONTINUE
985      RETURN
986  127 IDJ = 2-IDO
987      DO 130 J=2,IP
988         IDJ = IDJ+IDO
989         DO 129 K=1,L1
990            IDIJ = IDJ
991CDIR$ IVDEP
992            DO 128 I=4,IDO,2
993               IDIJ = IDIJ+2
994               C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
995               C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
996  128       CONTINUE
997  129    CONTINUE
998  130 CONTINUE
999      RETURN
1000      END
1001*
1002      SUBROUTINE PASSB2(IDO,L1,CC,CH,WA1)
1003CSE
1004      IMPLICIT REAL*8(A-H,O-Z)
1005C***BEGIN PROLOGUE  PASSB2
1006C***REFER TO  CFFTB
1007C***ROUTINES CALLED  (NONE)
1008C***END PROLOGUE  PASSB2
1009      DIMENSION       CC(IDO,2,L1)           ,CH(IDO,L1,2)           ,
1010     1                WA1(1)
1011C***FIRST EXECUTABLE STATEMENT  PASSB2
1012      IF (IDO .GT. 2) GO TO 102
1013      DO 101 K=1,L1
1014         CH(1,K,1) = CC(1,1,K)+CC(1,2,K)
1015         CH(1,K,2) = CC(1,1,K)-CC(1,2,K)
1016         CH(2,K,1) = CC(2,1,K)+CC(2,2,K)
1017         CH(2,K,2) = CC(2,1,K)-CC(2,2,K)
1018  101 CONTINUE
1019      RETURN
1020  102 IF(IDO/2.LT.L1) GO TO 105
1021      DO 104 K=1,L1
1022CDIR$ IVDEP
1023         DO 103 I=2,IDO,2
1024            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
1025            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
1026            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
1027            TI2 = CC(I,1,K)-CC(I,2,K)
1028            CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2
1029            CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2
1030  103    CONTINUE
1031  104 CONTINUE
1032      RETURN
1033  105 DO 107 I=2,IDO,2
1034CDIR$ IVDEP
1035         DO 106 K=1,L1
1036            CH(I-1,K,1) = CC(I-1,1,K)+CC(I-1,2,K)
1037            TR2 = CC(I-1,1,K)-CC(I-1,2,K)
1038            CH(I,K,1) = CC(I,1,K)+CC(I,2,K)
1039            TI2 = CC(I,1,K)-CC(I,2,K)
1040            CH(I,K,2) = WA1(I-1)*TI2+WA1(I)*TR2
1041            CH(I-1,K,2) = WA1(I-1)*TR2-WA1(I)*TI2
1042  106    CONTINUE
1043  107 CONTINUE
1044      RETURN
1045      END
1046*
1047      SUBROUTINE PASSB3(IDO,L1,CC,CH,WA1,WA2)
1048CSE
1049      IMPLICIT REAL*8(A-H,O-Z)
1050C***BEGIN PROLOGUE  PASSB3
1051C***REFER TO  CFFTB
1052C***ROUTINES CALLED  (NONE)
1053C***END PROLOGUE  PASSB3
1054      DIMENSION       CC(IDO,3,L1)           ,CH(IDO,L1,3)           ,
1055     1                WA1(1)     ,WA2(1)
1056      DATA TAUR,TAUI /-.5,.866025403784439/
1057C***FIRST EXECUTABLE STATEMENT  PASSB3
1058      IF (IDO .NE. 2) GO TO 102
1059      DO 101 K=1,L1
1060         TR2 = CC(1,2,K)+CC(1,3,K)
1061         CR2 = CC(1,1,K)+TAUR*TR2
1062         CH(1,K,1) = CC(1,1,K)+TR2
1063         TI2 = CC(2,2,K)+CC(2,3,K)
1064         CI2 = CC(2,1,K)+TAUR*TI2
1065         CH(2,K,1) = CC(2,1,K)+TI2
1066         CR3 = TAUI*(CC(1,2,K)-CC(1,3,K))
1067         CI3 = TAUI*(CC(2,2,K)-CC(2,3,K))
1068         CH(1,K,2) = CR2-CI3
1069         CH(1,K,3) = CR2+CI3
1070         CH(2,K,2) = CI2+CR3
1071         CH(2,K,3) = CI2-CR3
1072  101 CONTINUE
1073      RETURN
1074  102 IF(IDO/2.LT.L1) GO TO 105
1075      DO 104 K=1,L1
1076CDIR$ IVDEP
1077         DO 103 I=2,IDO,2
1078            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
1079            CR2 = CC(I-1,1,K)+TAUR*TR2
1080            CH(I-1,K,1) = CC(I-1,1,K)+TR2
1081            TI2 = CC(I,2,K)+CC(I,3,K)
1082            CI2 = CC(I,1,K)+TAUR*TI2
1083            CH(I,K,1) = CC(I,1,K)+TI2
1084            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
1085            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
1086            DR2 = CR2-CI3
1087            DR3 = CR2+CI3
1088            DI2 = CI2+CR3
1089            DI3 = CI2-CR3
1090            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
1091            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
1092            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
1093            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
1094  103    CONTINUE
1095  104 CONTINUE
1096      RETURN
1097  105 DO 107 I=2,IDO,2
1098CDIR$ IVDEP
1099         DO 106 K=1,L1
1100            TR2 = CC(I-1,2,K)+CC(I-1,3,K)
1101            CR2 = CC(I-1,1,K)+TAUR*TR2
1102            CH(I-1,K,1) = CC(I-1,1,K)+TR2
1103            TI2 = CC(I,2,K)+CC(I,3,K)
1104            CI2 = CC(I,1,K)+TAUR*TI2
1105            CH(I,K,1) = CC(I,1,K)+TI2
1106            CR3 = TAUI*(CC(I-1,2,K)-CC(I-1,3,K))
1107            CI3 = TAUI*(CC(I,2,K)-CC(I,3,K))
1108            DR2 = CR2-CI3
1109            DR3 = CR2+CI3
1110            DI2 = CI2+CR3
1111            DI3 = CI2-CR3
1112            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
1113            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
1114            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
1115            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
1116  106    CONTINUE
1117  107 CONTINUE
1118      RETURN
1119      END
1120*
1121      SUBROUTINE PASSB4(IDO,L1,CC,CH,WA1,WA2,WA3)
1122CSE
1123      IMPLICIT REAL*8(A-H,O-Z)
1124C***BEGIN PROLOGUE  PASSB4
1125C***REFER TO  CFFTB
1126C***ROUTINES CALLED  (NONE)
1127C***END PROLOGUE  PASSB4
1128      DIMENSION       CC(IDO,4,L1)           ,CH(IDO,L1,4)           ,
1129     1                WA1(1)     ,WA2(1)     ,WA3(1)
1130C***FIRST EXECUTABLE STATEMENT  PASSB4
1131      IF (IDO .NE. 2) GO TO 102
1132      DO 101 K=1,L1
1133         TI1 = CC(2,1,K)-CC(2,3,K)
1134         TI2 = CC(2,1,K)+CC(2,3,K)
1135         TR4 = CC(2,4,K)-CC(2,2,K)
1136         TI3 = CC(2,2,K)+CC(2,4,K)
1137         TR1 = CC(1,1,K)-CC(1,3,K)
1138         TR2 = CC(1,1,K)+CC(1,3,K)
1139         TI4 = CC(1,2,K)-CC(1,4,K)
1140         TR3 = CC(1,2,K)+CC(1,4,K)
1141         CH(1,K,1) = TR2+TR3
1142         CH(1,K,3) = TR2-TR3
1143         CH(2,K,1) = TI2+TI3
1144         CH(2,K,3) = TI2-TI3
1145         CH(1,K,2) = TR1+TR4
1146         CH(1,K,4) = TR1-TR4
1147         CH(2,K,2) = TI1+TI4
1148         CH(2,K,4) = TI1-TI4
1149  101 CONTINUE
1150      RETURN
1151  102 IF(IDO/2.LT.L1) GO TO 105
1152      DO 104 K=1,L1
1153CDIR$ IVDEP
1154         DO 103 I=2,IDO,2
1155            TI1 = CC(I,1,K)-CC(I,3,K)
1156            TI2 = CC(I,1,K)+CC(I,3,K)
1157            TI3 = CC(I,2,K)+CC(I,4,K)
1158            TR4 = CC(I,4,K)-CC(I,2,K)
1159            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
1160            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
1161            TI4 = CC(I-1,2,K)-CC(I-1,4,K)
1162            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
1163            CH(I-1,K,1) = TR2+TR3
1164            CR3 = TR2-TR3
1165            CH(I,K,1) = TI2+TI3
1166            CI3 = TI2-TI3
1167            CR2 = TR1+TR4
1168            CR4 = TR1-TR4
1169            CI2 = TI1+TI4
1170            CI4 = TI1-TI4
1171            CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2
1172            CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2
1173            CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3
1174            CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3
1175            CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4
1176            CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4
1177  103    CONTINUE
1178  104 CONTINUE
1179      RETURN
1180  105 DO 107 I=2,IDO,2
1181CDIR$ IVDEP
1182         DO 106 K=1,L1
1183            TI1 = CC(I,1,K)-CC(I,3,K)
1184            TI2 = CC(I,1,K)+CC(I,3,K)
1185            TI3 = CC(I,2,K)+CC(I,4,K)
1186            TR4 = CC(I,4,K)-CC(I,2,K)
1187            TR1 = CC(I-1,1,K)-CC(I-1,3,K)
1188            TR2 = CC(I-1,1,K)+CC(I-1,3,K)
1189            TI4 = CC(I-1,2,K)-CC(I-1,4,K)
1190            TR3 = CC(I-1,2,K)+CC(I-1,4,K)
1191            CH(I-1,K,1) = TR2+TR3
1192            CR3 = TR2-TR3
1193            CH(I,K,1) = TI2+TI3
1194            CI3 = TI2-TI3
1195            CR2 = TR1+TR4
1196            CR4 = TR1-TR4
1197            CI2 = TI1+TI4
1198            CI4 = TI1-TI4
1199            CH(I-1,K,2) = WA1(I-1)*CR2-WA1(I)*CI2
1200            CH(I,K,2) = WA1(I-1)*CI2+WA1(I)*CR2
1201            CH(I-1,K,3) = WA2(I-1)*CR3-WA2(I)*CI3
1202            CH(I,K,3) = WA2(I-1)*CI3+WA2(I)*CR3
1203            CH(I-1,K,4) = WA3(I-1)*CR4-WA3(I)*CI4
1204            CH(I,K,4) = WA3(I-1)*CI4+WA3(I)*CR4
1205  106    CONTINUE
1206  107 CONTINUE
1207      RETURN
1208      END
1209*
1210      SUBROUTINE PASSB5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1211CSE
1212      IMPLICIT REAL*8(A-H,O-Z)
1213C***BEGIN PROLOGUE  PASSB5
1214C***REFER TO  CFFTB
1215C***ROUTINES CALLED  (NONE)
1216C***END PROLOGUE  PASSB5
1217      DIMENSION       CC(IDO,5,L1)           ,CH(IDO,L1,5)           ,
1218     1                WA1(1)     ,WA2(1)     ,WA3(1)     ,WA4(1)
1219      DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
1220     1-.809016994374947,.587785252292473/
1221C***FIRST EXECUTABLE STATEMENT  PASSB5
1222      IF (IDO .NE. 2) GO TO 102
1223      DO 101 K=1,L1
1224         TI5 = CC(2,2,K)-CC(2,5,K)
1225         TI2 = CC(2,2,K)+CC(2,5,K)
1226         TI4 = CC(2,3,K)-CC(2,4,K)
1227         TI3 = CC(2,3,K)+CC(2,4,K)
1228         TR5 = CC(1,2,K)-CC(1,5,K)
1229         TR2 = CC(1,2,K)+CC(1,5,K)
1230         TR4 = CC(1,3,K)-CC(1,4,K)
1231         TR3 = CC(1,3,K)+CC(1,4,K)
1232         CH(1,K,1) = CC(1,1,K)+TR2+TR3
1233         CH(2,K,1) = CC(2,1,K)+TI2+TI3
1234         CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
1235         CI2 = CC(2,1,K)+TR11*TI2+TR12*TI3
1236         CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
1237         CI3 = CC(2,1,K)+TR12*TI2+TR11*TI3
1238         CR5 = TI11*TR5+TI12*TR4
1239         CI5 = TI11*TI5+TI12*TI4
1240         CR4 = TI12*TR5-TI11*TR4
1241         CI4 = TI12*TI5-TI11*TI4
1242         CH(1,K,2) = CR2-CI5
1243         CH(1,K,5) = CR2+CI5
1244         CH(2,K,2) = CI2+CR5
1245         CH(2,K,3) = CI3+CR4
1246         CH(1,K,3) = CR3-CI4
1247         CH(1,K,4) = CR3+CI4
1248         CH(2,K,4) = CI3-CR4
1249         CH(2,K,5) = CI2-CR5
1250  101 CONTINUE
1251      RETURN
1252  102 IF(IDO/2.LT.L1) GO TO 105
1253      DO 104 K=1,L1
1254CDIR$ IVDEP
1255         DO 103 I=2,IDO,2
1256            TI5 = CC(I,2,K)-CC(I,5,K)
1257            TI2 = CC(I,2,K)+CC(I,5,K)
1258            TI4 = CC(I,3,K)-CC(I,4,K)
1259            TI3 = CC(I,3,K)+CC(I,4,K)
1260            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
1261            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
1262            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
1263            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
1264            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
1265            CH(I,K,1) = CC(I,1,K)+TI2+TI3
1266            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
1267            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
1268            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
1269            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
1270            CR5 = TI11*TR5+TI12*TR4
1271            CI5 = TI11*TI5+TI12*TI4
1272            CR4 = TI12*TR5-TI11*TR4
1273            CI4 = TI12*TI5-TI11*TI4
1274            DR3 = CR3-CI4
1275            DR4 = CR3+CI4
1276            DI3 = CI3+CR4
1277            DI4 = CI3-CR4
1278            DR5 = CR2+CI5
1279            DR2 = CR2-CI5
1280            DI5 = CI2-CR5
1281            DI2 = CI2+CR5
1282            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
1283            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
1284            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
1285            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
1286            CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4
1287            CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4
1288            CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5
1289            CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5
1290  103    CONTINUE
1291  104 CONTINUE
1292      RETURN
1293  105 DO 107 I=2,IDO,2
1294CDIR$ IVDEP
1295         DO 106 K=1,L1
1296            TI5 = CC(I,2,K)-CC(I,5,K)
1297            TI2 = CC(I,2,K)+CC(I,5,K)
1298            TI4 = CC(I,3,K)-CC(I,4,K)
1299            TI3 = CC(I,3,K)+CC(I,4,K)
1300            TR5 = CC(I-1,2,K)-CC(I-1,5,K)
1301            TR2 = CC(I-1,2,K)+CC(I-1,5,K)
1302            TR4 = CC(I-1,3,K)-CC(I-1,4,K)
1303            TR3 = CC(I-1,3,K)+CC(I-1,4,K)
1304            CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
1305            CH(I,K,1) = CC(I,1,K)+TI2+TI3
1306            CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
1307            CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
1308            CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
1309            CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
1310            CR5 = TI11*TR5+TI12*TR4
1311            CI5 = TI11*TI5+TI12*TI4
1312            CR4 = TI12*TR5-TI11*TR4
1313            CI4 = TI12*TI5-TI11*TI4
1314            DR3 = CR3-CI4
1315            DR4 = CR3+CI4
1316            DI3 = CI3+CR4
1317            DI4 = CI3-CR4
1318            DR5 = CR2+CI5
1319            DR2 = CR2-CI5
1320            DI5 = CI2-CR5
1321            DI2 = CI2+CR5
1322            CH(I-1,K,2) = WA1(I-1)*DR2-WA1(I)*DI2
1323            CH(I,K,2) = WA1(I-1)*DI2+WA1(I)*DR2
1324            CH(I-1,K,3) = WA2(I-1)*DR3-WA2(I)*DI3
1325            CH(I,K,3) = WA2(I-1)*DI3+WA2(I)*DR3
1326            CH(I-1,K,4) = WA3(I-1)*DR4-WA3(I)*DI4
1327            CH(I,K,4) = WA3(I-1)*DI4+WA3(I)*DR4
1328            CH(I-1,K,5) = WA4(I-1)*DR5-WA4(I)*DI5
1329            CH(I,K,5) = WA4(I-1)*DI5+WA4(I)*DR5
1330  106    CONTINUE
1331  107 CONTINUE
1332      RETURN
1333      END
1334