/[MITgcm]/MITgcm/pkg/zonal_filt/fftpack.F
ViewVC logotype

Annotation of /MITgcm/pkg/zonal_filt/fftpack.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph

Revision 1.3 - (hide annotations) (download)
Sun Feb 4 14:38:50 2001 UTC (23 years, 3 months ago) by cnh
Branch: MAIN
CVS Tags: checkpoint38, c37_adj, checkpoint39, checkpoint37, checkpoint36, checkpoint35
Branch point for: pre38
Changes since 1.2: +2 -0 lines 
Made sure each .F and .h file had
the CVS keywords Header and Name at its start.
Most had header but very few currently have Name, so
lots of changes!
1 cnh 1.3 C $Header: $
2     C $Name: $
3 adcroft 1.2 C SUBROUTINE RFFTB (N,R,WSAVE)
4     C DIMENSION R(1) ,WSAVE(1)
5     C IF (N .EQ. 1) RETURN
6     C CALL RFFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
7     C RETURN
8     C END
9     SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1)
10     DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
11     1 WA1(1)
12     DO 101 K=1,L1
13     CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
14     CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
15     101 CONTINUE
16     IF (IDO-2) 107,105,102
17     102 IDP2 = IDO+2
18     DO 104 K=1,L1
19     DO 103 I=3,IDO,2
20     IC = IDP2-I
21     CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
22     TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
23     CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
24     TI2 = CC(I,1,K)+CC(IC,2,K)
25     CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
26     CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
27     103 CONTINUE
28     104 CONTINUE
29     IF (MOD(IDO,2) .EQ. 1) RETURN
30     105 DO 106 K=1,L1
31     CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
32     CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
33     106 CONTINUE
34     107 RETURN
35     END
36     SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2)
37     DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
38     1 WA1(1) ,WA2(1)
39     DATA TAUR,TAUI /-.5,.866025403784439/
40     DO 101 K=1,L1
41     TR2 = CC(IDO,2,K)+CC(IDO,2,K)
42     CR2 = CC(1,1,K)+TAUR*TR2
43     CH(1,K,1) = CC(1,1,K)+TR2
44     CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
45     CH(1,K,2) = CR2-CI3
46     CH(1,K,3) = CR2+CI3
47     101 CONTINUE
48     IF (IDO .EQ. 1) RETURN
49     IDP2 = IDO+2
50     DO 103 K=1,L1
51     DO 102 I=3,IDO,2
52     IC = IDP2-I
53     TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
54     CR2 = CC(I-1,1,K)+TAUR*TR2
55     CH(I-1,K,1) = CC(I-1,1,K)+TR2
56     TI2 = CC(I,3,K)-CC(IC,2,K)
57     CI2 = CC(I,1,K)+TAUR*TI2
58     CH(I,K,1) = CC(I,1,K)+TI2
59     CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
60     CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
61     DR2 = CR2-CI3
62     DR3 = CR2+CI3
63     DI2 = CI2+CR3
64     DI3 = CI2-CR3
65     CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
66     CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
67     CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
68     CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
69     102 CONTINUE
70     103 CONTINUE
71     RETURN
72     END
73     SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
74     DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
75     1 WA1(1) ,WA2(1) ,WA3(1)
76     DATA SQRT2 /1.414213562373095/
77     DO 101 K=1,L1
78     TR1 = CC(1,1,K)-CC(IDO,4,K)
79     TR2 = CC(1,1,K)+CC(IDO,4,K)
80     TR3 = CC(IDO,2,K)+CC(IDO,2,K)
81     TR4 = CC(1,3,K)+CC(1,3,K)
82     CH(1,K,1) = TR2+TR3
83     CH(1,K,2) = TR1-TR4
84     CH(1,K,3) = TR2-TR3
85     CH(1,K,4) = TR1+TR4
86     101 CONTINUE
87     IF (IDO-2) 107,105,102
88     102 IDP2 = IDO+2
89     DO 104 K=1,L1
90     DO 103 I=3,IDO,2
91     IC = IDP2-I
92     TI1 = CC(I,1,K)+CC(IC,4,K)
93     TI2 = CC(I,1,K)-CC(IC,4,K)
94     TI3 = CC(I,3,K)-CC(IC,2,K)
95     TR4 = CC(I,3,K)+CC(IC,2,K)
96     TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
97     TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
98     TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
99     TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
100     CH(I-1,K,1) = TR2+TR3
101     CR3 = TR2-TR3
102     CH(I,K,1) = TI2+TI3
103     CI3 = TI2-TI3
104     CR2 = TR1-TR4
105     CR4 = TR1+TR4
106     CI2 = TI1+TI4
107     CI4 = TI1-TI4
108     CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
109     CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
110     CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
111     CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
112     CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
113     CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
114     103 CONTINUE
115     104 CONTINUE
116     IF (MOD(IDO,2) .EQ. 1) RETURN
117     105 CONTINUE
118     DO 106 K=1,L1
119     TI1 = CC(1,2,K)+CC(1,4,K)
120     TI2 = CC(1,4,K)-CC(1,2,K)
121     TR1 = CC(IDO,1,K)-CC(IDO,3,K)
122     TR2 = CC(IDO,1,K)+CC(IDO,3,K)
123     CH(IDO,K,1) = TR2+TR2
124     CH(IDO,K,2) = SQRT2*(TR1-TI1)
125     CH(IDO,K,3) = TI2+TI2
126     CH(IDO,K,4) = -SQRT2*(TR1+TI1)
127     106 CONTINUE
128     107 RETURN
129     END
130     SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
131     DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
132     1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
133     DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
134     1-.809016994374947,.587785252292473/
135     DO 101 K=1,L1
136     TI5 = CC(1,3,K)+CC(1,3,K)
137     TI4 = CC(1,5,K)+CC(1,5,K)
138     TR2 = CC(IDO,2,K)+CC(IDO,2,K)
139     TR3 = CC(IDO,4,K)+CC(IDO,4,K)
140     CH(1,K,1) = CC(1,1,K)+TR2+TR3
141     CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
142     CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
143     CI5 = TI11*TI5+TI12*TI4
144     CI4 = TI12*TI5-TI11*TI4
145     CH(1,K,2) = CR2-CI5
146     CH(1,K,3) = CR3-CI4
147     CH(1,K,4) = CR3+CI4
148     CH(1,K,5) = CR2+CI5
149     101 CONTINUE
150     IF (IDO .EQ. 1) RETURN
151     IDP2 = IDO+2
152     DO 103 K=1,L1
153     DO 102 I=3,IDO,2
154     IC = IDP2-I
155     TI5 = CC(I,3,K)+CC(IC,2,K)
156     TI2 = CC(I,3,K)-CC(IC,2,K)
157     TI4 = CC(I,5,K)+CC(IC,4,K)
158     TI3 = CC(I,5,K)-CC(IC,4,K)
159     TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
160     TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
161     TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
162     TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
163     CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
164     CH(I,K,1) = CC(I,1,K)+TI2+TI3
165     CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
166     CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
167     CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
168     CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
169     CR5 = TI11*TR5+TI12*TR4
170     CI5 = TI11*TI5+TI12*TI4
171     CR4 = TI12*TR5-TI11*TR4
172     CI4 = TI12*TI5-TI11*TI4
173     DR3 = CR3-CI4
174     DR4 = CR3+CI4
175     DI3 = CI3+CR4
176     DI4 = CI3-CR4
177     DR5 = CR2+CI5
178     DR2 = CR2-CI5
179     DI5 = CI2-CR5
180     DI2 = CI2+CR5
181     CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
182     CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
183     CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
184     CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
185     CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
186     CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
187     CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
188     CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
189     102 CONTINUE
190     103 CONTINUE
191     RETURN
192     END
193     SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
194     DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
195     1 C1(IDO,L1,IP) ,C2(IDL1,IP),
196     2 CH2(IDL1,IP) ,WA(1)
197     DATA TPI/6.28318530717959/
198     ARG = TPI/FLOAT(IP)
199     DCP = COS(ARG)
200     DSP = SIN(ARG)
201     IDP2 = IDO+2
202     NBD = (IDO-1)/2
203     IPP2 = IP+2
204     IPPH = (IP+1)/2
205     IF (IDO .LT. L1) GO TO 103
206     DO 102 K=1,L1
207     DO 101 I=1,IDO
208     CH(I,K,1) = CC(I,1,K)
209     101 CONTINUE
210     102 CONTINUE
211     GO TO 106
212     103 DO 105 I=1,IDO
213     DO 104 K=1,L1
214     CH(I,K,1) = CC(I,1,K)
215     104 CONTINUE
216     105 CONTINUE
217     106 DO 108 J=2,IPPH
218     JC = IPP2-J
219     J2 = J+J
220     DO 107 K=1,L1
221     CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
222     CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
223     107 CONTINUE
224     108 CONTINUE
225     IF (IDO .EQ. 1) GO TO 116
226     IF (NBD .LT. L1) GO TO 112
227     DO 111 J=2,IPPH
228     JC = IPP2-J
229     DO 110 K=1,L1
230     DO 109 I=3,IDO,2
231     IC = IDP2-I
232     CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
233     CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
234     CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
235     CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
236     109 CONTINUE
237     110 CONTINUE
238     111 CONTINUE
239     GO TO 116
240     112 DO 115 J=2,IPPH
241     JC = IPP2-J
242     DO 114 I=3,IDO,2
243     IC = IDP2-I
244     DO 113 K=1,L1
245     CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
246     CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
247     CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
248     CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
249     113 CONTINUE
250     114 CONTINUE
251     115 CONTINUE
252     116 AR1 = 1.
253     AI1 = 0.
254     DO 120 L=2,IPPH
255     LC = IPP2-L
256     AR1H = DCP*AR1-DSP*AI1
257     AI1 = DCP*AI1+DSP*AR1
258     AR1 = AR1H
259     DO 117 IK=1,IDL1
260     C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
261     C2(IK,LC) = AI1*CH2(IK,IP)
262     117 CONTINUE
263     DC2 = AR1
264     DS2 = AI1
265     AR2 = AR1
266     AI2 = AI1
267     DO 119 J=3,IPPH
268     JC = IPP2-J
269     AR2H = DC2*AR2-DS2*AI2
270     AI2 = DC2*AI2+DS2*AR2
271     AR2 = AR2H
272     DO 118 IK=1,IDL1
273     C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
274     C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
275     118 CONTINUE
276     119 CONTINUE
277     120 CONTINUE
278     DO 122 J=2,IPPH
279     DO 121 IK=1,IDL1
280     CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
281     121 CONTINUE
282     122 CONTINUE
283     DO 124 J=2,IPPH
284     JC = IPP2-J
285     DO 123 K=1,L1
286     CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
287     CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
288     123 CONTINUE
289     124 CONTINUE
290     IF (IDO .EQ. 1) GO TO 132
291     IF (NBD .LT. L1) GO TO 128
292     DO 127 J=2,IPPH
293     JC = IPP2-J
294     DO 126 K=1,L1
295     DO 125 I=3,IDO,2
296     CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
297     CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
298     CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
299     CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
300     125 CONTINUE
301     126 CONTINUE
302     127 CONTINUE
303     GO TO 132
304     128 DO 131 J=2,IPPH
305     JC = IPP2-J
306     DO 130 I=3,IDO,2
307     DO 129 K=1,L1
308     CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
309     CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
310     CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
311     CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
312     129 CONTINUE
313     130 CONTINUE
314     131 CONTINUE
315     132 CONTINUE
316     IF (IDO .EQ. 1) RETURN
317     DO 133 IK=1,IDL1
318     C2(IK,1) = CH2(IK,1)
319     133 CONTINUE
320     DO 135 J=2,IP
321     DO 134 K=1,L1
322     C1(1,K,J) = CH(1,K,J)
323     134 CONTINUE
324     135 CONTINUE
325     IF (NBD .GT. L1) GO TO 139
326     IS = -IDO
327     DO 138 J=2,IP
328     IS = IS+IDO
329     IDIJ = IS
330     DO 137 I=3,IDO,2
331     IDIJ = IDIJ+2
332     DO 136 K=1,L1
333     C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
334     C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
335     136 CONTINUE
336     137 CONTINUE
337     138 CONTINUE
338     GO TO 143
339     139 IS = -IDO
340     DO 142 J=2,IP
341     IS = IS+IDO
342     DO 141 K=1,L1
343     IDIJ = IS
344     DO 140 I=3,IDO,2
345     IDIJ = IDIJ+2
346     C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
347     C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
348     140 CONTINUE
349     141 CONTINUE
350     142 CONTINUE
351     143 RETURN
352     END
353     SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC)
354     DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
355     NF = IFAC(2)
356     NA = 0
357     L1 = 1
358     IW = 1
359     DO 116 K1=1,NF
360     IP = IFAC(K1+2)
361     L2 = IP*L1
362     IDO = N/L2
363     IDL1 = IDO*L1
364     IF (IP .NE. 4) GO TO 103
365     IX2 = IW+IDO
366     IX3 = IX2+IDO
367     IF (NA .NE. 0) GO TO 101
368     CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
369     GO TO 102
370     101 CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
371     102 NA = 1-NA
372     GO TO 115
373     103 IF (IP .NE. 2) GO TO 106
374     IF (NA .NE. 0) GO TO 104
375     CALL RADB2 (IDO,L1,C,CH,WA(IW))
376     GO TO 105
377     104 CALL RADB2 (IDO,L1,CH,C,WA(IW))
378     105 NA = 1-NA
379     GO TO 115
380     106 IF (IP .NE. 3) GO TO 109
381     IX2 = IW+IDO
382     IF (NA .NE. 0) GO TO 107
383     CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
384     GO TO 108
385     107 CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
386     108 NA = 1-NA
387     GO TO 115
388     109 IF (IP .NE. 5) GO TO 112
389     IX2 = IW+IDO
390     IX3 = IX2+IDO
391     IX4 = IX3+IDO
392     IF (NA .NE. 0) GO TO 110
393     CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
394     GO TO 111
395     110 CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
396     111 NA = 1-NA
397     GO TO 115
398     112 IF (NA .NE. 0) GO TO 113
399     CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
400     GO TO 114
401     113 CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
402     114 IF (IDO .EQ. 1) NA = 1-NA
403     115 L1 = L2
404     IW = IW+(IP-1)*IDO
405     116 CONTINUE
406     IF (NA .EQ. 0) RETURN
407     DO 117 I=1,N
408     C(I) = CH(I)
409     117 CONTINUE
410     RETURN
411     END
412     C SUBROUTINE RFFTF (N,R,WSAVE)
413     C DIMENSION R(1) ,WSAVE(1)
414     C IF (N .EQ. 1) RETURN
415     C CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
416     C RETURN
417     C END
418     SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1)
419     DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
420     1 WA1(1)
421     DO 101 K=1,L1
422     CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
423     CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
424     101 CONTINUE
425     IF (IDO-2) 107,105,102
426     102 IDP2 = IDO+2
427     DO 104 K=1,L1
428     DO 103 I=3,IDO,2
429     IC = IDP2-I
430     TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
431     TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
432     CH(I,1,K) = CC(I,K,1)+TI2
433     CH(IC,2,K) = TI2-CC(I,K,1)
434     CH(I-1,1,K) = CC(I-1,K,1)+TR2
435     CH(IC-1,2,K) = CC(I-1,K,1)-TR2
436     103 CONTINUE
437     104 CONTINUE
438     IF (MOD(IDO,2) .EQ. 1) RETURN
439     105 DO 106 K=1,L1
440     CH(1,2,K) = -CC(IDO,K,2)
441     CH(IDO,1,K) = CC(IDO,K,1)
442     106 CONTINUE
443     107 RETURN
444     END
445     SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)
446     DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
447     1 WA1(1) ,WA2(1)
448     DATA TAUR,TAUI /-.5,.866025403784439/
449     DO 101 K=1,L1
450     CR2 = CC(1,K,2)+CC(1,K,3)
451     CH(1,1,K) = CC(1,K,1)+CR2
452     CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
453     CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
454     101 CONTINUE
455     IF (IDO .EQ. 1) RETURN
456     IDP2 = IDO+2
457     DO 103 K=1,L1
458     DO 102 I=3,IDO,2
459     IC = IDP2-I
460     DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
461     DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
462     DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
463     DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
464     CR2 = DR2+DR3
465     CI2 = DI2+DI3
466     CH(I-1,1,K) = CC(I-1,K,1)+CR2
467     CH(I,1,K) = CC(I,K,1)+CI2
468     TR2 = CC(I-1,K,1)+TAUR*CR2
469     TI2 = CC(I,K,1)+TAUR*CI2
470     TR3 = TAUI*(DI2-DI3)
471     TI3 = TAUI*(DR3-DR2)
472     CH(I-1,3,K) = TR2+TR3
473     CH(IC-1,2,K) = TR2-TR3
474     CH(I,3,K) = TI2+TI3
475     CH(IC,2,K) = TI3-TI2
476     102 CONTINUE
477     103 CONTINUE
478     RETURN
479     END
480     SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
481     DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
482     1 WA1(1) ,WA2(1) ,WA3(1)
483     DATA HSQT2 /.7071067811865475/
484     DO 101 K=1,L1
485     TR1 = CC(1,K,2)+CC(1,K,4)
486     TR2 = CC(1,K,1)+CC(1,K,3)
487     CH(1,1,K) = TR1+TR2
488     CH(IDO,4,K) = TR2-TR1
489     CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
490     CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
491     101 CONTINUE
492     IF (IDO-2) 107,105,102
493     102 IDP2 = IDO+2
494     DO 104 K=1,L1
495     DO 103 I=3,IDO,2
496     IC = IDP2-I
497     CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
498     CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
499     CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
500     CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
501     CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
502     CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
503     TR1 = CR2+CR4
504     TR4 = CR4-CR2
505     TI1 = CI2+CI4
506     TI4 = CI2-CI4
507     TI2 = CC(I,K,1)+CI3
508     TI3 = CC(I,K,1)-CI3
509     TR2 = CC(I-1,K,1)+CR3
510     TR3 = CC(I-1,K,1)-CR3
511     CH(I-1,1,K) = TR1+TR2
512     CH(IC-1,4,K) = TR2-TR1
513     CH(I,1,K) = TI1+TI2
514     CH(IC,4,K) = TI1-TI2
515     CH(I-1,3,K) = TI4+TR3
516     CH(IC-1,2,K) = TR3-TI4
517     CH(I,3,K) = TR4+TI3
518     CH(IC,2,K) = TR4-TI3
519     103 CONTINUE
520     104 CONTINUE
521     IF (MOD(IDO,2) .EQ. 1) RETURN
522     105 CONTINUE
523     DO 106 K=1,L1
524     TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
525     TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
526     CH(IDO,1,K) = TR1+CC(IDO,K,1)
527     CH(IDO,3,K) = CC(IDO,K,1)-TR1
528     CH(1,2,K) = TI1-CC(IDO,K,3)
529     CH(1,4,K) = TI1+CC(IDO,K,3)
530     106 CONTINUE
531     107 RETURN
532     END
533     SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
534     DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
535     1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
536     DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
537     1-.809016994374947,.587785252292473/
538     DO 101 K=1,L1
539     CR2 = CC(1,K,5)+CC(1,K,2)
540     CI5 = CC(1,K,5)-CC(1,K,2)
541     CR3 = CC(1,K,4)+CC(1,K,3)
542     CI4 = CC(1,K,4)-CC(1,K,3)
543     CH(1,1,K) = CC(1,K,1)+CR2+CR3
544     CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
545     CH(1,3,K) = TI11*CI5+TI12*CI4
546     CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
547     CH(1,5,K) = TI12*CI5-TI11*CI4
548     101 CONTINUE
549     IF (IDO .EQ. 1) RETURN
550     IDP2 = IDO+2
551     DO 103 K=1,L1
552     DO 102 I=3,IDO,2
553     IC = IDP2-I
554     DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
555     DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
556     DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
557     DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
558     DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
559     DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
560     DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
561     DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
562     CR2 = DR2+DR5
563     CI5 = DR5-DR2
564     CR5 = DI2-DI5
565     CI2 = DI2+DI5
566     CR3 = DR3+DR4
567     CI4 = DR4-DR3
568     CR4 = DI3-DI4
569     CI3 = DI3+DI4
570     CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
571     CH(I,1,K) = CC(I,K,1)+CI2+CI3
572     TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
573     TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
574     TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
575     TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
576     TR5 = TI11*CR5+TI12*CR4
577     TI5 = TI11*CI5+TI12*CI4
578     TR4 = TI12*CR5-TI11*CR4
579     TI4 = TI12*CI5-TI11*CI4
580     CH(I-1,3,K) = TR2+TR5
581     CH(IC-1,2,K) = TR2-TR5
582     CH(I,3,K) = TI2+TI5
583     CH(IC,2,K) = TI5-TI2
584     CH(I-1,5,K) = TR3+TR4
585     CH(IC-1,4,K) = TR3-TR4
586     CH(I,5,K) = TI3+TI4
587     CH(IC,4,K) = TI4-TI3
588     102 CONTINUE
589     103 CONTINUE
590     RETURN
591     END
592     SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
593     DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
594     1 C1(IDO,L1,IP) ,C2(IDL1,IP),
595     2 CH2(IDL1,IP) ,WA(1)
596     DATA TPI/6.28318530717959/
597     ARG = TPI/FLOAT(IP)
598     DCP = COS(ARG)
599     DSP = SIN(ARG)
600     IPPH = (IP+1)/2
601     IPP2 = IP+2
602     IDP2 = IDO+2
603     NBD = (IDO-1)/2
604     IF (IDO .EQ. 1) GO TO 119
605     DO 101 IK=1,IDL1
606     CH2(IK,1) = C2(IK,1)
607     101 CONTINUE
608     DO 103 J=2,IP
609     DO 102 K=1,L1
610     CH(1,K,J) = C1(1,K,J)
611     102 CONTINUE
612     103 CONTINUE
613     IF (NBD .GT. L1) GO TO 107
614     IS = -IDO
615     DO 106 J=2,IP
616     IS = IS+IDO
617     IDIJ = IS
618     DO 105 I=3,IDO,2
619     IDIJ = IDIJ+2
620     DO 104 K=1,L1
621     CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
622     CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
623     104 CONTINUE
624     105 CONTINUE
625     106 CONTINUE
626     GO TO 111
627     107 IS = -IDO
628     DO 110 J=2,IP
629     IS = IS+IDO
630     DO 109 K=1,L1
631     IDIJ = IS
632     DO 108 I=3,IDO,2
633     IDIJ = IDIJ+2
634     CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
635     CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
636     108 CONTINUE
637     109 CONTINUE
638     110 CONTINUE
639     111 IF (NBD .LT. L1) GO TO 115
640     DO 114 J=2,IPPH
641     JC = IPP2-J
642     DO 113 K=1,L1
643     DO 112 I=3,IDO,2
644     C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
645     C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
646     C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
647     C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
648     112 CONTINUE
649     113 CONTINUE
650     114 CONTINUE
651     GO TO 121
652     115 DO 118 J=2,IPPH
653     JC = IPP2-J
654     DO 117 I=3,IDO,2
655     DO 116 K=1,L1
656     C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
657     C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
658     C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
659     C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
660     116 CONTINUE
661     117 CONTINUE
662     118 CONTINUE
663     GO TO 121
664     119 DO 120 IK=1,IDL1
665     C2(IK,1) = CH2(IK,1)
666     120 CONTINUE
667     121 DO 123 J=2,IPPH
668     JC = IPP2-J
669     DO 122 K=1,L1
670     C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
671     C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
672     122 CONTINUE
673     123 CONTINUE
674     C
675     AR1 = 1.
676     AI1 = 0.
677     DO 127 L=2,IPPH
678     LC = IPP2-L
679     AR1H = DCP*AR1-DSP*AI1
680     AI1 = DCP*AI1+DSP*AR1
681     AR1 = AR1H
682     DO 124 IK=1,IDL1
683     CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
684     CH2(IK,LC) = AI1*C2(IK,IP)
685     124 CONTINUE
686     DC2 = AR1
687     DS2 = AI1
688     AR2 = AR1
689     AI2 = AI1
690     DO 126 J=3,IPPH
691     JC = IPP2-J
692     AR2H = DC2*AR2-DS2*AI2
693     AI2 = DC2*AI2+DS2*AR2
694     AR2 = AR2H
695     DO 125 IK=1,IDL1
696     CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
697     CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
698     125 CONTINUE
699     126 CONTINUE
700     127 CONTINUE
701     DO 129 J=2,IPPH
702     DO 128 IK=1,IDL1
703     CH2(IK,1) = CH2(IK,1)+C2(IK,J)
704     128 CONTINUE
705     129 CONTINUE
706     C
707     IF (IDO .LT. L1) GO TO 132
708     DO 131 K=1,L1
709     DO 130 I=1,IDO
710     CC(I,1,K) = CH(I,K,1)
711     130 CONTINUE
712     131 CONTINUE
713     GO TO 135
714     132 DO 134 I=1,IDO
715     DO 133 K=1,L1
716     CC(I,1,K) = CH(I,K,1)
717     133 CONTINUE
718     134 CONTINUE
719     135 DO 137 J=2,IPPH
720     JC = IPP2-J
721     J2 = J+J
722     DO 136 K=1,L1
723     CC(IDO,J2-2,K) = CH(1,K,J)
724     CC(1,J2-1,K) = CH(1,K,JC)
725     136 CONTINUE
726     137 CONTINUE
727     IF (IDO .EQ. 1) RETURN
728     IF (NBD .LT. L1) GO TO 141
729     DO 140 J=2,IPPH
730     JC = IPP2-J
731     J2 = J+J
732     DO 139 K=1,L1
733     DO 138 I=3,IDO,2
734     IC = IDP2-I
735     CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
736     CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
737     CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
738     CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
739     138 CONTINUE
740     139 CONTINUE
741     140 CONTINUE
742     RETURN
743     141 DO 144 J=2,IPPH
744     JC = IPP2-J
745     J2 = J+J
746     DO 143 I=3,IDO,2
747     IC = IDP2-I
748     DO 142 K=1,L1
749     CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
750     CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
751     CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
752     CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
753     142 CONTINUE
754     143 CONTINUE
755     144 CONTINUE
756     RETURN
757     END
758     SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC)
759     DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
760     NF = IFAC(2)
761     NA = 1
762     L2 = N
763     IW = N
764     DO 111 K1=1,NF
765     KH = NF-K1
766     IP = IFAC(KH+3)
767     L1 = L2/IP
768     IDO = N/L2
769     IDL1 = IDO*L1
770     IW = IW-(IP-1)*IDO
771     NA = 1-NA
772     IF (IP .NE. 4) GO TO 102
773     IX2 = IW+IDO
774     IX3 = IX2+IDO
775     IF (NA .NE. 0) GO TO 101
776     CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
777     GO TO 110
778     101 CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
779     GO TO 110
780     102 IF (IP .NE. 2) GO TO 104
781     IF (NA .NE. 0) GO TO 103
782     CALL RADF2 (IDO,L1,C,CH,WA(IW))
783     GO TO 110
784     103 CALL RADF2 (IDO,L1,CH,C,WA(IW))
785     GO TO 110
786     104 IF (IP .NE. 3) GO TO 106
787     IX2 = IW+IDO
788     IF (NA .NE. 0) GO TO 105
789     CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
790     GO TO 110
791     105 CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
792     GO TO 110
793     106 IF (IP .NE. 5) GO TO 108
794     IX2 = IW+IDO
795     IX3 = IX2+IDO
796     IX4 = IX3+IDO
797     IF (NA .NE. 0) GO TO 107
798     CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
799     GO TO 110
800     107 CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
801     GO TO 110
802     108 IF (IDO .EQ. 1) NA = 1-NA
803     IF (NA .NE. 0) GO TO 109
804     CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
805     NA = 1
806     GO TO 110
807     109 CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
808     NA = 0
809     110 L2 = L1
810     111 CONTINUE
811     IF (NA .EQ. 1) RETURN
812     DO 112 I=1,N
813     C(I) = CH(I)
814     112 CONTINUE
815     RETURN
816     END
817     C SUBROUTINE RFFTI (N,WSAVE)
818     C DIMENSION WSAVE(1)
819     C IF (N .EQ. 1) RETURN
820     C CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))
821     C RETURN
822     C END
823     SUBROUTINE RFFTI1 (N,WA,IFAC)
824     DIMENSION WA(1) ,IFAC(1) ,NTRYH(4)
825     DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
826     NL = N
827     NF = 0
828     J = 0
829     101 J = J+1
830     IF (J-4) 102,102,103
831     102 NTRY = NTRYH(J)
832     GO TO 104
833     103 NTRY = NTRY+2
834     104 NQ = NL/NTRY
835     NR = NL-NTRY*NQ
836     IF (NR) 101,105,101
837     105 NF = NF+1
838     IFAC(NF+2) = NTRY
839     NL = NQ
840     IF (NTRY .NE. 2) GO TO 107
841     IF (NF .EQ. 1) GO TO 107
842     DO 106 I=2,NF
843     IB = NF-I+2
844     IFAC(IB+2) = IFAC(IB+1)
845     106 CONTINUE
846     IFAC(3) = 2
847     107 IF (NL .NE. 1) GO TO 104
848     IFAC(1) = N
849     IFAC(2) = NF
850     TPI = 6.28318530717959
851     ARGH = TPI/FLOAT(N)
852     IS = 0
853     NFM1 = NF-1
854     L1 = 1
855     IF (NFM1 .EQ. 0) RETURN
856     DO 110 K1=1,NFM1
857     IP = IFAC(K1+2)
858     LD = 0
859     L2 = L1*IP
860     IDO = N/L2
861     IPM = IP-1
862     DO 109 J=1,IPM
863     LD = LD+L1
864     I = IS
865     ARGLD = FLOAT(LD)*ARGH
866     FI = 0.
867     DO 108 II=3,IDO,2
868     I = I+2
869     FI = FI+1.
870     ARG = FI*ARGLD
871     WA(I-1) = COS(ARG)
872     WA(I) = SIN(ARG)
873     108 CONTINUE
874     IS = IS+IDO
875     109 CONTINUE
876     L1 = L2
877     110 CONTINUE
878     RETURN
879     END
880     C SUBROUTINE R8FFTB (N,R,WSAVE)
881     C implicit real*8 (A-H,O-Z)
882     C DIMENSION R(1) ,WSAVE(1)
883     C IF (N .EQ. 1) RETURN
884     C CALL R8FFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
885     C RETURN
886     C END
887     SUBROUTINE R8ADB2 (IDO,L1,CC,CH,WA1)
888     implicit real*8 (A-H,O-Z)
889     DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
890     1 WA1(1)
891     DO 101 K=1,L1
892     CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
893     CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
894     101 CONTINUE
895     IF (IDO-2) 107,105,102
896     102 IDP2 = IDO+2
897     DO 104 K=1,L1
898     DO 103 I=3,IDO,2
899     IC = IDP2-I
900     CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
901     TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
902     CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
903     TI2 = CC(I,1,K)+CC(IC,2,K)
904     CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
905     CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
906     103 CONTINUE
907     104 CONTINUE
908     IF (MOD(IDO,2) .EQ. 1) RETURN
909     105 DO 106 K=1,L1
910     CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
911     CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
912     106 CONTINUE
913     107 RETURN
914     END
915     SUBROUTINE R8ADB3 (IDO,L1,CC,CH,WA1,WA2)
916     implicit real*8 (A-H,O-Z)
917     DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
918     1 WA1(1) ,WA2(1)
919     DATA TAUR,TAUI /-.5,.866025403784439/
920     DO 101 K=1,L1
921     TR2 = CC(IDO,2,K)+CC(IDO,2,K)
922     CR2 = CC(1,1,K)+TAUR*TR2
923     CH(1,K,1) = CC(1,1,K)+TR2
924     CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
925     CH(1,K,2) = CR2-CI3
926     CH(1,K,3) = CR2+CI3
927     101 CONTINUE
928     IF (IDO .EQ. 1) RETURN
929     IDP2 = IDO+2
930     DO 103 K=1,L1
931     DO 102 I=3,IDO,2
932     IC = IDP2-I
933     TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
934     CR2 = CC(I-1,1,K)+TAUR*TR2
935     CH(I-1,K,1) = CC(I-1,1,K)+TR2
936     TI2 = CC(I,3,K)-CC(IC,2,K)
937     CI2 = CC(I,1,K)+TAUR*TI2
938     CH(I,K,1) = CC(I,1,K)+TI2
939     CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
940     CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
941     DR2 = CR2-CI3
942     DR3 = CR2+CI3
943     DI2 = CI2+CR3
944     DI3 = CI2-CR3
945     CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
946     CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
947     CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
948     CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
949     102 CONTINUE
950     103 CONTINUE
951     RETURN
952     END
953     SUBROUTINE R8ADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
954     implicit real*8 (A-H,O-Z)
955     DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
956     1 WA1(1) ,WA2(1) ,WA3(1)
957     DATA SQRT2 /1.414213562373095/
958     DO 101 K=1,L1
959     TR1 = CC(1,1,K)-CC(IDO,4,K)
960     TR2 = CC(1,1,K)+CC(IDO,4,K)
961     TR3 = CC(IDO,2,K)+CC(IDO,2,K)
962     TR4 = CC(1,3,K)+CC(1,3,K)
963     CH(1,K,1) = TR2+TR3
964     CH(1,K,2) = TR1-TR4
965     CH(1,K,3) = TR2-TR3
966     CH(1,K,4) = TR1+TR4
967     101 CONTINUE
968     IF (IDO-2) 107,105,102
969     102 IDP2 = IDO+2
970     DO 104 K=1,L1
971     DO 103 I=3,IDO,2
972     IC = IDP2-I
973     TI1 = CC(I,1,K)+CC(IC,4,K)
974     TI2 = CC(I,1,K)-CC(IC,4,K)
975     TI3 = CC(I,3,K)-CC(IC,2,K)
976     TR4 = CC(I,3,K)+CC(IC,2,K)
977     TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
978     TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
979     TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
980     TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
981     CH(I-1,K,1) = TR2+TR3
982     CR3 = TR2-TR3
983     CH(I,K,1) = TI2+TI3
984     CI3 = TI2-TI3
985     CR2 = TR1-TR4
986     CR4 = TR1+TR4
987     CI2 = TI1+TI4
988     CI4 = TI1-TI4
989     CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
990     CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
991     CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
992     CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
993     CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
994     CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
995     103 CONTINUE
996     104 CONTINUE
997     IF (MOD(IDO,2) .EQ. 1) RETURN
998     105 CONTINUE
999     DO 106 K=1,L1
1000     TI1 = CC(1,2,K)+CC(1,4,K)
1001     TI2 = CC(1,4,K)-CC(1,2,K)
1002     TR1 = CC(IDO,1,K)-CC(IDO,3,K)
1003     TR2 = CC(IDO,1,K)+CC(IDO,3,K)
1004     CH(IDO,K,1) = TR2+TR2
1005     CH(IDO,K,2) = SQRT2*(TR1-TI1)
1006     CH(IDO,K,3) = TI2+TI2
1007     CH(IDO,K,4) = -SQRT2*(TR1+TI1)
1008     106 CONTINUE
1009     107 RETURN
1010     END
1011     SUBROUTINE R8ADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1012     implicit real*8 (A-H,O-Z)
1013     DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
1014     1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
1015     DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
1016     1-.809016994374947,.587785252292473/
1017     DO 101 K=1,L1
1018     TI5 = CC(1,3,K)+CC(1,3,K)
1019     TI4 = CC(1,5,K)+CC(1,5,K)
1020     TR2 = CC(IDO,2,K)+CC(IDO,2,K)
1021     TR3 = CC(IDO,4,K)+CC(IDO,4,K)
1022     CH(1,K,1) = CC(1,1,K)+TR2+TR3
1023     CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
1024     CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
1025     CI5 = TI11*TI5+TI12*TI4
1026     CI4 = TI12*TI5-TI11*TI4
1027     CH(1,K,2) = CR2-CI5
1028     CH(1,K,3) = CR3-CI4
1029     CH(1,K,4) = CR3+CI4
1030     CH(1,K,5) = CR2+CI5
1031     101 CONTINUE
1032     IF (IDO .EQ. 1) RETURN
1033     IDP2 = IDO+2
1034     DO 103 K=1,L1
1035     DO 102 I=3,IDO,2
1036     IC = IDP2-I
1037     TI5 = CC(I,3,K)+CC(IC,2,K)
1038     TI2 = CC(I,3,K)-CC(IC,2,K)
1039     TI4 = CC(I,5,K)+CC(IC,4,K)
1040     TI3 = CC(I,5,K)-CC(IC,4,K)
1041     TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
1042     TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
1043     TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
1044     TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
1045     CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
1046     CH(I,K,1) = CC(I,1,K)+TI2+TI3
1047     CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
1048     CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
1049     CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
1050     CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
1051     CR5 = TI11*TR5+TI12*TR4
1052     CI5 = TI11*TI5+TI12*TI4
1053     CR4 = TI12*TR5-TI11*TR4
1054     CI4 = TI12*TI5-TI11*TI4
1055     DR3 = CR3-CI4
1056     DR4 = CR3+CI4
1057     DI3 = CI3+CR4
1058     DI4 = CI3-CR4
1059     DR5 = CR2+CI5
1060     DR2 = CR2-CI5
1061     DI5 = CI2-CR5
1062     DI2 = CI2+CR5
1063     CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
1064     CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
1065     CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
1066     CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
1067     CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
1068     CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
1069     CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
1070     CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
1071     102 CONTINUE
1072     103 CONTINUE
1073     RETURN
1074     END
1075     SUBROUTINE R8ADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1076     implicit real*8 (A-H,O-Z)
1077     DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1078     1 C1(IDO,L1,IP) ,C2(IDL1,IP),
1079     2 CH2(IDL1,IP) ,WA(1)
1080     DATA TPI/6.28318530717959/
1081     ARG = TPI/FLOAT(IP)
1082     DCP = COS(ARG)
1083     DSP = SIN(ARG)
1084     IDP2 = IDO+2
1085     NBD = (IDO-1)/2
1086     IPP2 = IP+2
1087     IPPH = (IP+1)/2
1088     IF (IDO .LT. L1) GO TO 103
1089     DO 102 K=1,L1
1090     DO 101 I=1,IDO
1091     CH(I,K,1) = CC(I,1,K)
1092     101 CONTINUE
1093     102 CONTINUE
1094     GO TO 106
1095     103 DO 105 I=1,IDO
1096     DO 104 K=1,L1
1097     CH(I,K,1) = CC(I,1,K)
1098     104 CONTINUE
1099     105 CONTINUE
1100     106 DO 108 J=2,IPPH
1101     JC = IPP2-J
1102     J2 = J+J
1103     DO 107 K=1,L1
1104     CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
1105     CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
1106     107 CONTINUE
1107     108 CONTINUE
1108     IF (IDO .EQ. 1) GO TO 116
1109     IF (NBD .LT. L1) GO TO 112
1110     DO 111 J=2,IPPH
1111     JC = IPP2-J
1112     DO 110 K=1,L1
1113     DO 109 I=3,IDO,2
1114     IC = IDP2-I
1115     CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
1116     CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
1117     CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
1118     CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
1119     109 CONTINUE
1120     110 CONTINUE
1121     111 CONTINUE
1122     GO TO 116
1123     112 DO 115 J=2,IPPH
1124     JC = IPP2-J
1125     DO 114 I=3,IDO,2
1126     IC = IDP2-I
1127     DO 113 K=1,L1
1128     CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
1129     CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
1130     CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
1131     CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
1132     113 CONTINUE
1133     114 CONTINUE
1134     115 CONTINUE
1135     116 AR1 = 1.
1136     AI1 = 0.
1137     DO 120 L=2,IPPH
1138     LC = IPP2-L
1139     AR1H = DCP*AR1-DSP*AI1
1140     AI1 = DCP*AI1+DSP*AR1
1141     AR1 = AR1H
1142     DO 117 IK=1,IDL1
1143     C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
1144     C2(IK,LC) = AI1*CH2(IK,IP)
1145     117 CONTINUE
1146     DC2 = AR1
1147     DS2 = AI1
1148     AR2 = AR1
1149     AI2 = AI1
1150     DO 119 J=3,IPPH
1151     JC = IPP2-J
1152     AR2H = DC2*AR2-DS2*AI2
1153     AI2 = DC2*AI2+DS2*AR2
1154     AR2 = AR2H
1155     DO 118 IK=1,IDL1
1156     C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
1157     C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
1158     118 CONTINUE
1159     119 CONTINUE
1160     120 CONTINUE
1161     DO 122 J=2,IPPH
1162     DO 121 IK=1,IDL1
1163     CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
1164     121 CONTINUE
1165     122 CONTINUE
1166     DO 124 J=2,IPPH
1167     JC = IPP2-J
1168     DO 123 K=1,L1
1169     CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
1170     CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
1171     123 CONTINUE
1172     124 CONTINUE
1173     IF (IDO .EQ. 1) GO TO 132
1174     IF (NBD .LT. L1) GO TO 128
1175     DO 127 J=2,IPPH
1176     JC = IPP2-J
1177     DO 126 K=1,L1
1178     DO 125 I=3,IDO,2
1179     CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
1180     CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
1181     CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
1182     CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
1183     125 CONTINUE
1184     126 CONTINUE
1185     127 CONTINUE
1186     GO TO 132
1187     128 DO 131 J=2,IPPH
1188     JC = IPP2-J
1189     DO 130 I=3,IDO,2
1190     DO 129 K=1,L1
1191     CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
1192     CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
1193     CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
1194     CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
1195     129 CONTINUE
1196     130 CONTINUE
1197     131 CONTINUE
1198     132 CONTINUE
1199     IF (IDO .EQ. 1) RETURN
1200     DO 133 IK=1,IDL1
1201     C2(IK,1) = CH2(IK,1)
1202     133 CONTINUE
1203     DO 135 J=2,IP
1204     DO 134 K=1,L1
1205     C1(1,K,J) = CH(1,K,J)
1206     134 CONTINUE
1207     135 CONTINUE
1208     IF (NBD .GT. L1) GO TO 139
1209     IS = -IDO
1210     DO 138 J=2,IP
1211     IS = IS+IDO
1212     IDIJ = IS
1213     DO 137 I=3,IDO,2
1214     IDIJ = IDIJ+2
1215     DO 136 K=1,L1
1216     C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
1217     C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
1218     136 CONTINUE
1219     137 CONTINUE
1220     138 CONTINUE
1221     GO TO 143
1222     139 IS = -IDO
1223     DO 142 J=2,IP
1224     IS = IS+IDO
1225     DO 141 K=1,L1
1226     IDIJ = IS
1227     DO 140 I=3,IDO,2
1228     IDIJ = IDIJ+2
1229     C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
1230     C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
1231     140 CONTINUE
1232     141 CONTINUE
1233     142 CONTINUE
1234     143 RETURN
1235     END
1236     SUBROUTINE R8FFTB1 (N,C,CH,WA,IFAC)
1237     implicit real*8 (A-H,O-Z)
1238     DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
1239     NF = IFAC(2)
1240     NA = 0
1241     L1 = 1
1242     IW = 1
1243     DO 116 K1=1,NF
1244     IP = IFAC(K1+2)
1245     L2 = IP*L1
1246     IDO = N/L2
1247     IDL1 = IDO*L1
1248     IF (IP .NE. 4) GO TO 103
1249     IX2 = IW+IDO
1250     IX3 = IX2+IDO
1251     IF (NA .NE. 0) GO TO 101
1252     CALL R8ADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
1253     GO TO 102
1254     101 CALL R8ADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
1255     102 NA = 1-NA
1256     GO TO 115
1257     103 IF (IP .NE. 2) GO TO 106
1258     IF (NA .NE. 0) GO TO 104
1259     CALL R8ADB2 (IDO,L1,C,CH,WA(IW))
1260     GO TO 105
1261     104 CALL R8ADB2 (IDO,L1,CH,C,WA(IW))
1262     105 NA = 1-NA
1263     GO TO 115
1264     106 IF (IP .NE. 3) GO TO 109
1265     IX2 = IW+IDO
1266     IF (NA .NE. 0) GO TO 107
1267     CALL R8ADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
1268     GO TO 108
1269     107 CALL R8ADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
1270     108 NA = 1-NA
1271     GO TO 115
1272     109 IF (IP .NE. 5) GO TO 112
1273     IX2 = IW+IDO
1274     IX3 = IX2+IDO
1275     IX4 = IX3+IDO
1276     IF (NA .NE. 0) GO TO 110
1277     CALL R8ADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1278     GO TO 111
1279     110 CALL R8ADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1280     111 NA = 1-NA
1281     GO TO 115
1282     112 IF (NA .NE. 0) GO TO 113
1283     CALL R8ADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
1284     GO TO 114
1285     113 CALL R8ADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
1286     114 IF (IDO .EQ. 1) NA = 1-NA
1287     115 L1 = L2
1288     IW = IW+(IP-1)*IDO
1289     116 CONTINUE
1290     IF (NA .EQ. 0) RETURN
1291     DO 117 I=1,N
1292     C(I) = CH(I)
1293     117 CONTINUE
1294     RETURN
1295     END
1296     C SUBROUTINE R8FFTF (N,R,WSAVE)
1297     C implicit real*8 (A-H,O-Z)
1298     C DIMENSION R(1) ,WSAVE(1)
1299     C IF (N .EQ. 1) RETURN
1300     C CALL R8FFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
1301     C RETURN
1302     C END
1303     SUBROUTINE R8ADF2 (IDO,L1,CC,CH,WA1)
1304     implicit real*8 (A-H,O-Z)
1305     DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
1306     1 WA1(1)
1307     DO 101 K=1,L1
1308     CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
1309     CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
1310     101 CONTINUE
1311     IF (IDO-2) 107,105,102
1312     102 IDP2 = IDO+2
1313     DO 104 K=1,L1
1314     DO 103 I=3,IDO,2
1315     IC = IDP2-I
1316     TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1317     TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1318     CH(I,1,K) = CC(I,K,1)+TI2
1319     CH(IC,2,K) = TI2-CC(I,K,1)
1320     CH(I-1,1,K) = CC(I-1,K,1)+TR2
1321     CH(IC-1,2,K) = CC(I-1,K,1)-TR2
1322     103 CONTINUE
1323     104 CONTINUE
1324     IF (MOD(IDO,2) .EQ. 1) RETURN
1325     105 DO 106 K=1,L1
1326     CH(1,2,K) = -CC(IDO,K,2)
1327     CH(IDO,1,K) = CC(IDO,K,1)
1328     106 CONTINUE
1329     107 RETURN
1330     END
1331     SUBROUTINE R8ADF3 (IDO,L1,CC,CH,WA1,WA2)
1332     implicit real*8 (A-H,O-Z)
1333     DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
1334     1 WA1(1) ,WA2(1)
1335     DATA TAUR,TAUI /-.5,.866025403784439/
1336     DO 101 K=1,L1
1337     CR2 = CC(1,K,2)+CC(1,K,3)
1338     CH(1,1,K) = CC(1,K,1)+CR2
1339     CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
1340     CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
1341     101 CONTINUE
1342     IF (IDO .EQ. 1) RETURN
1343     IDP2 = IDO+2
1344     DO 103 K=1,L1
1345     DO 102 I=3,IDO,2
1346     IC = IDP2-I
1347     DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1348     DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1349     DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1350     DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1351     CR2 = DR2+DR3
1352     CI2 = DI2+DI3
1353     CH(I-1,1,K) = CC(I-1,K,1)+CR2
1354     CH(I,1,K) = CC(I,K,1)+CI2
1355     TR2 = CC(I-1,K,1)+TAUR*CR2
1356     TI2 = CC(I,K,1)+TAUR*CI2
1357     TR3 = TAUI*(DI2-DI3)
1358     TI3 = TAUI*(DR3-DR2)
1359     CH(I-1,3,K) = TR2+TR3
1360     CH(IC-1,2,K) = TR2-TR3
1361     CH(I,3,K) = TI2+TI3
1362     CH(IC,2,K) = TI3-TI2
1363     102 CONTINUE
1364     103 CONTINUE
1365     RETURN
1366     END
1367     SUBROUTINE R8ADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1368     implicit real*8 (A-H,O-Z)
1369     DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
1370     1 WA1(1) ,WA2(1) ,WA3(1)
1371     DATA HSQT2 /.7071067811865475/
1372     DO 101 K=1,L1
1373     TR1 = CC(1,K,2)+CC(1,K,4)
1374     TR2 = CC(1,K,1)+CC(1,K,3)
1375     CH(1,1,K) = TR1+TR2
1376     CH(IDO,4,K) = TR2-TR1
1377     CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
1378     CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
1379     101 CONTINUE
1380     IF (IDO-2) 107,105,102
1381     102 IDP2 = IDO+2
1382     DO 104 K=1,L1
1383     DO 103 I=3,IDO,2
1384     IC = IDP2-I
1385     CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1386     CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1387     CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1388     CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1389     CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
1390     CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
1391     TR1 = CR2+CR4
1392     TR4 = CR4-CR2
1393     TI1 = CI2+CI4
1394     TI4 = CI2-CI4
1395     TI2 = CC(I,K,1)+CI3
1396     TI3 = CC(I,K,1)-CI3
1397     TR2 = CC(I-1,K,1)+CR3
1398     TR3 = CC(I-1,K,1)-CR3
1399     CH(I-1,1,K) = TR1+TR2
1400     CH(IC-1,4,K) = TR2-TR1
1401     CH(I,1,K) = TI1+TI2
1402     CH(IC,4,K) = TI1-TI2
1403     CH(I-1,3,K) = TI4+TR3
1404     CH(IC-1,2,K) = TR3-TI4
1405     CH(I,3,K) = TR4+TI3
1406     CH(IC,2,K) = TR4-TI3
1407     103 CONTINUE
1408     104 CONTINUE
1409     IF (MOD(IDO,2) .EQ. 1) RETURN
1410     105 CONTINUE
1411     DO 106 K=1,L1
1412     TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
1413     TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
1414     CH(IDO,1,K) = TR1+CC(IDO,K,1)
1415     CH(IDO,3,K) = CC(IDO,K,1)-TR1
1416     CH(1,2,K) = TI1-CC(IDO,K,3)
1417     CH(1,4,K) = TI1+CC(IDO,K,3)
1418     106 CONTINUE
1419     107 RETURN
1420     END
1421     SUBROUTINE R8ADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1422     implicit real*8 (A-H,O-Z)
1423     DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
1424     1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
1425     DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
1426     1-.809016994374947,.587785252292473/
1427     DO 101 K=1,L1
1428     CR2 = CC(1,K,5)+CC(1,K,2)
1429     CI5 = CC(1,K,5)-CC(1,K,2)
1430     CR3 = CC(1,K,4)+CC(1,K,3)
1431     CI4 = CC(1,K,4)-CC(1,K,3)
1432     CH(1,1,K) = CC(1,K,1)+CR2+CR3
1433     CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
1434     CH(1,3,K) = TI11*CI5+TI12*CI4
1435     CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
1436     CH(1,5,K) = TI12*CI5-TI11*CI4
1437     101 CONTINUE
1438     IF (IDO .EQ. 1) RETURN
1439     IDP2 = IDO+2
1440     DO 103 K=1,L1
1441     DO 102 I=3,IDO,2
1442     IC = IDP2-I
1443     DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
1444     DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
1445     DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
1446     DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
1447     DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
1448     DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
1449     DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
1450     DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
1451     CR2 = DR2+DR5
1452     CI5 = DR5-DR2
1453     CR5 = DI2-DI5
1454     CI2 = DI2+DI5
1455     CR3 = DR3+DR4
1456     CI4 = DR4-DR3
1457     CR4 = DI3-DI4
1458     CI3 = DI3+DI4
1459     CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
1460     CH(I,1,K) = CC(I,K,1)+CI2+CI3
1461     TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
1462     TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
1463     TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
1464     TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
1465     TR5 = TI11*CR5+TI12*CR4
1466     TI5 = TI11*CI5+TI12*CI4
1467     TR4 = TI12*CR5-TI11*CR4
1468     TI4 = TI12*CI5-TI11*CI4
1469     CH(I-1,3,K) = TR2+TR5
1470     CH(IC-1,2,K) = TR2-TR5
1471     CH(I,3,K) = TI2+TI5
1472     CH(IC,2,K) = TI5-TI2
1473     CH(I-1,5,K) = TR3+TR4
1474     CH(IC-1,4,K) = TR3-TR4
1475     CH(I,5,K) = TI3+TI4
1476     CH(IC,4,K) = TI4-TI3
1477     102 CONTINUE
1478     103 CONTINUE
1479     RETURN
1480     END
1481     SUBROUTINE R8ADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1482     implicit real*8 (A-H,O-Z)
1483     DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
1484     1 C1(IDO,L1,IP) ,C2(IDL1,IP),
1485     2 CH2(IDL1,IP) ,WA(1)
1486     DATA TPI/6.28318530717959/
1487     ARG = TPI/FLOAT(IP)
1488     DCP = COS(ARG)
1489     DSP = SIN(ARG)
1490     IPPH = (IP+1)/2
1491     IPP2 = IP+2
1492     IDP2 = IDO+2
1493     NBD = (IDO-1)/2
1494     IF (IDO .EQ. 1) GO TO 119
1495     DO 101 IK=1,IDL1
1496     CH2(IK,1) = C2(IK,1)
1497     101 CONTINUE
1498     DO 103 J=2,IP
1499     DO 102 K=1,L1
1500     CH(1,K,J) = C1(1,K,J)
1501     102 CONTINUE
1502     103 CONTINUE
1503     IF (NBD .GT. L1) GO TO 107
1504     IS = -IDO
1505     DO 106 J=2,IP
1506     IS = IS+IDO
1507     IDIJ = IS
1508     DO 105 I=3,IDO,2
1509     IDIJ = IDIJ+2
1510     DO 104 K=1,L1
1511     CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
1512     CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
1513     104 CONTINUE
1514     105 CONTINUE
1515     106 CONTINUE
1516     GO TO 111
1517     107 IS = -IDO
1518     DO 110 J=2,IP
1519     IS = IS+IDO
1520     DO 109 K=1,L1
1521     IDIJ = IS
1522     DO 108 I=3,IDO,2
1523     IDIJ = IDIJ+2
1524     CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
1525     CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
1526     108 CONTINUE
1527     109 CONTINUE
1528     110 CONTINUE
1529     111 IF (NBD .LT. L1) GO TO 115
1530     DO 114 J=2,IPPH
1531     JC = IPP2-J
1532     DO 113 K=1,L1
1533     DO 112 I=3,IDO,2
1534     C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
1535     C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
1536     C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
1537     C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
1538     112 CONTINUE
1539     113 CONTINUE
1540     114 CONTINUE
1541     GO TO 121
1542     115 DO 118 J=2,IPPH
1543     JC = IPP2-J
1544     DO 117 I=3,IDO,2
1545     DO 116 K=1,L1
1546     C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
1547     C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
1548     C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
1549     C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
1550     116 CONTINUE
1551     117 CONTINUE
1552     118 CONTINUE
1553     GO TO 121
1554     119 DO 120 IK=1,IDL1
1555     C2(IK,1) = CH2(IK,1)
1556     120 CONTINUE
1557     121 DO 123 J=2,IPPH
1558     JC = IPP2-J
1559     DO 122 K=1,L1
1560     C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
1561     C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
1562     122 CONTINUE
1563     123 CONTINUE
1564     C
1565     AR1 = 1.
1566     AI1 = 0.
1567     DO 127 L=2,IPPH
1568     LC = IPP2-L
1569     AR1H = DCP*AR1-DSP*AI1
1570     AI1 = DCP*AI1+DSP*AR1
1571     AR1 = AR1H
1572     DO 124 IK=1,IDL1
1573     CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
1574     CH2(IK,LC) = AI1*C2(IK,IP)
1575     124 CONTINUE
1576     DC2 = AR1
1577     DS2 = AI1
1578     AR2 = AR1
1579     AI2 = AI1
1580     DO 126 J=3,IPPH
1581     JC = IPP2-J
1582     AR2H = DC2*AR2-DS2*AI2
1583     AI2 = DC2*AI2+DS2*AR2
1584     AR2 = AR2H
1585     DO 125 IK=1,IDL1
1586     CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
1587     CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
1588     125 CONTINUE
1589     126 CONTINUE
1590     127 CONTINUE
1591     DO 129 J=2,IPPH
1592     DO 128 IK=1,IDL1
1593     CH2(IK,1) = CH2(IK,1)+C2(IK,J)
1594     128 CONTINUE
1595     129 CONTINUE
1596     C
1597     IF (IDO .LT. L1) GO TO 132
1598     DO 131 K=1,L1
1599     DO 130 I=1,IDO
1600     CC(I,1,K) = CH(I,K,1)
1601     130 CONTINUE
1602     131 CONTINUE
1603     GO TO 135
1604     132 DO 134 I=1,IDO
1605     DO 133 K=1,L1
1606     CC(I,1,K) = CH(I,K,1)
1607     133 CONTINUE
1608     134 CONTINUE
1609     135 DO 137 J=2,IPPH
1610     JC = IPP2-J
1611     J2 = J+J
1612     DO 136 K=1,L1
1613     CC(IDO,J2-2,K) = CH(1,K,J)
1614     CC(1,J2-1,K) = CH(1,K,JC)
1615     136 CONTINUE
1616     137 CONTINUE
1617     IF (IDO .EQ. 1) RETURN
1618     IF (NBD .LT. L1) GO TO 141
1619     DO 140 J=2,IPPH
1620     JC = IPP2-J
1621     J2 = J+J
1622     DO 139 K=1,L1
1623     DO 138 I=3,IDO,2
1624     IC = IDP2-I
1625     CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1626     CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1627     CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1628     CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1629     138 CONTINUE
1630     139 CONTINUE
1631     140 CONTINUE
1632     RETURN
1633     141 DO 144 J=2,IPPH
1634     JC = IPP2-J
1635     J2 = J+J
1636     DO 143 I=3,IDO,2
1637     IC = IDP2-I
1638     DO 142 K=1,L1
1639     CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1640     CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1641     CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1642     CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1643     142 CONTINUE
1644     143 CONTINUE
1645     144 CONTINUE
1646     RETURN
1647     END
1648     SUBROUTINE R8FFTF1 (N,C,CH,WA,IFAC)
1649     implicit real*8 (A-H,O-Z)
1650     DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(1)
1651     NF = IFAC(2)
1652     NA = 1
1653     L2 = N
1654     IW = N
1655     DO 111 K1=1,NF
1656     KH = NF-K1
1657     IP = IFAC(KH+3)
1658     L1 = L2/IP
1659     IDO = N/L2
1660     IDL1 = IDO*L1
1661     IW = IW-(IP-1)*IDO
1662     NA = 1-NA
1663     IF (IP .NE. 4) GO TO 102
1664     IX2 = IW+IDO
1665     IX3 = IX2+IDO
1666     IF (NA .NE. 0) GO TO 101
1667     CALL R8ADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
1668     GO TO 110
1669     101 CALL R8ADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
1670     GO TO 110
1671     102 IF (IP .NE. 2) GO TO 104
1672     IF (NA .NE. 0) GO TO 103
1673     CALL R8ADF2 (IDO,L1,C,CH,WA(IW))
1674     GO TO 110
1675     103 CALL R8ADF2 (IDO,L1,CH,C,WA(IW))
1676     GO TO 110
1677     104 IF (IP .NE. 3) GO TO 106
1678     IX2 = IW+IDO
1679     IF (NA .NE. 0) GO TO 105
1680     CALL R8ADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
1681     GO TO 110
1682     105 CALL R8ADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
1683     GO TO 110
1684     106 IF (IP .NE. 5) GO TO 108
1685     IX2 = IW+IDO
1686     IX3 = IX2+IDO
1687     IX4 = IX3+IDO
1688     IF (NA .NE. 0) GO TO 107
1689     CALL R8ADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1690     GO TO 110
1691     107 CALL R8ADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
1692     GO TO 110
1693     108 IF (IDO .EQ. 1) NA = 1-NA
1694     IF (NA .NE. 0) GO TO 109
1695     CALL R8ADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
1696     NA = 1
1697     GO TO 110
1698     109 CALL R8ADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
1699     NA = 0
1700     110 L2 = L1
1701     111 CONTINUE
1702     IF (NA .EQ. 1) RETURN
1703     DO 112 I=1,N
1704     C(I) = CH(I)
1705     112 CONTINUE
1706     RETURN
1707     END
1708     C SUBROUTINE R8FFTI (N,WSAVE)
1709     C implicit real*8 (A-H,O-Z)
1710     C DIMENSION WSAVE(1)
1711     C IF (N .EQ. 1) RETURN
1712     C CALL R8FFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))
1713     C RETURN
1714     C END
1715     SUBROUTINE R8FFTI1 (N,WA,IFAC)
1716     implicit real*8 (A-H,O-Z)
1717     DIMENSION WA(1) ,IFAC(1) ,NTRYH(4)
1718     DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
1719     NL = N
1720     NF = 0
1721     J = 0
1722     101 J = J+1
1723     IF (J-4) 102,102,103
1724     102 NTRY = NTRYH(J)
1725     GO TO 104
1726     103 NTRY = NTRY+2
1727     104 NQ = NL/NTRY
1728     NR = NL-NTRY*NQ
1729     IF (NR) 101,105,101
1730     105 NF = NF+1
1731     IFAC(NF+2) = NTRY
1732     NL = NQ
1733     IF (NTRY .NE. 2) GO TO 107
1734     IF (NF .EQ. 1) GO TO 107
1735     DO 106 I=2,NF
1736     IB = NF-I+2
1737     IFAC(IB+2) = IFAC(IB+1)
1738     106 CONTINUE
1739     IFAC(3) = 2
1740     107 IF (NL .NE. 1) GO TO 104
1741     IFAC(1) = N
1742     IFAC(2) = NF
1743     TPI = 6.28318530717959
1744     ARGH = TPI/FLOAT(N)
1745     IS = 0
1746     NFM1 = NF-1
1747     L1 = 1
1748     IF (NFM1 .EQ. 0) RETURN
1749     DO 110 K1=1,NFM1
1750     IP = IFAC(K1+2)
1751     LD = 0
1752     L2 = L1*IP
1753     IDO = N/L2
1754     IPM = IP-1
1755     DO 109 J=1,IPM
1756     LD = LD+L1
1757     I = IS
1758     ARGLD = FLOAT(LD)*ARGH
1759     FI = 0.
1760     DO 108 II=3,IDO,2
1761     I = I+2
1762     FI = FI+1.
1763     ARG = FI*ARGLD
1764     WA(I-1) = COS(ARG)
1765     WA(I) = SIN(ARG)
1766     108 CONTINUE
1767     IS = IS+IDO
1768     109 CONTINUE
1769     L1 = L2
1770     110 CONTINUE
1771     RETURN
1772     END
  ViewVC Help
Powered by ViewVC 1.1.22