File: C:\NOAA\NEMS_11731\src\atmos\gfs\dyn\fftpack.f
1
2 SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale,
3 & table,n1,wrk,n2,z,nz)
4
5 implicit none
6 integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
7 real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
8
9 IF (init.ne.0) THEN
10 CALL rffti(n,table)
11 ELSE
12
13 DO j=1,m
14 y(1,j)=x(1,j)
15 DO i=2,n
16 y(i,j)=x(i+1,j)
17 ENDDO
18 CALL rfftb(n,y(1,j),table)
19 DO i=1,n
20 y(i,j)=scale*y(i,j)
21 ENDDO
22 ENDDO
23 ENDIF
24
25 RETURN
26 END
27
28
29
30 SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys)
31
32 implicit none
33 integer isign,n,isys,i
34 real scale,x(*),y(*),table(*),work(*)
35
36 IF (isign.eq.0) THEN
37 CALL rffti(n,table)
38 ENDIF
39 IF (isign.eq.1) THEN
40 y(1)=x(1)
41 DO i=2,n
42 y(i)=x(i+1)
43 ENDDO
44 CALL rfftb(n,y,table)
45 DO i=1,n
46 y(i)=scale*y(i)
47 ENDDO
48 ENDIF
49
50 RETURN
51 END
52
53
54
55 SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale,
56 & table,n1,wrk,n2,z,nz)
57
58 implicit none
59 integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
60 real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
61
62 IF (init.ne.0) THEN
63 CALL rffti(n,table)
64 ELSE
65 DO j=1,m
66 DO i=1,n
67 y(i,j)=x(i,j)
68 ENDDO
69 CALL rfftf(n,y(1,j),table)
70 DO i=1,n
71 y(i,j)=scale*y(i,j)
72 ENDDO
73 DO i=n,2,-1
74 y(i+1,j)=y(i,j)
75 ENDDO
76 y(2,j)=0.
77 ENDDO
78 ENDIF
79
80 RETURN
81 END
82
83
84
85 SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys)
86
87 implicit none
88 integer isign,n,isys,i
89 real scale,x(*),y(*),table(*),work(*)
90
91 IF (isign.eq.0) THEN
92 CALL rffti(n,table)
93 ENDIF
94 IF (isign.eq.-1) THEN
95 DO i=1,n
96 y(i)=x(i)
97 ENDDO
98 CALL rfftf(n,y,table)
99 DO i=1,n
100 y(i)=scale*y(i)
101 ENDDO
102 DO i=n,2,-1
103 y(i+1)=y(i)
104 ENDDO
105 y(2)=0.
106 ENDIF
107
108 RETURN
109 END
110
111
112
113
114
115
116
117
118
119 SUBROUTINE RFFTF (N,R,WSAVE)
120 DIMENSION R(1) ,WSAVE(1)
121 IF (N .EQ. 1) RETURN
122 CALL RFFTF1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
123 RETURN
124 END
125 SUBROUTINE RFFTB (N,R,WSAVE)
126 DIMENSION R(1) ,WSAVE(1)
127 IF (N .EQ. 1) RETURN
128 CALL RFFTB1 (N,R,WSAVE,WSAVE(N+1),WSAVE(2*N+1))
129 RETURN
130 END
131 SUBROUTINE RFFTI (N,WSAVE)
132 DIMENSION WSAVE(1)
133 IF (N .EQ. 1) RETURN
134 CALL RFFTI1 (N,WSAVE(N+1),WSAVE(2*N+1))
135 RETURN
136 END
137 SUBROUTINE RFFTB1 (N,C,CH,WA,IFAC)
138 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(*)
139 NF = IFAC(2)
140 NA = 0
141 L1 = 1
142 IW = 1
143 DO 116 K1=1,NF
144 IP = IFAC(K1+2)
145 L2 = IP*L1
146 IDO = N/L2
147 IDL1 = IDO*L1
148 IF (IP .NE. 4) GO TO 103
149 IX2 = IW+IDO
150 IX3 = IX2+IDO
151 IF (NA .NE. 0) GO TO 101
152 CALL RADB4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
153 GO TO 102
154 101 CALL RADB4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
155 102 NA = 1-NA
156 GO TO 115
157 103 IF (IP .NE. 2) GO TO 106
158 IF (NA .NE. 0) GO TO 104
159 CALL RADB2 (IDO,L1,C,CH,WA(IW))
160 GO TO 105
161 104 CALL RADB2 (IDO,L1,CH,C,WA(IW))
162 105 NA = 1-NA
163 GO TO 115
164 106 IF (IP .NE. 3) GO TO 109
165 IX2 = IW+IDO
166 IF (NA .NE. 0) GO TO 107
167 CALL RADB3 (IDO,L1,C,CH,WA(IW),WA(IX2))
168 GO TO 108
169 107 CALL RADB3 (IDO,L1,CH,C,WA(IW),WA(IX2))
170 108 NA = 1-NA
171 GO TO 115
172 109 IF (IP .NE. 5) GO TO 112
173 IX2 = IW+IDO
174 IX3 = IX2+IDO
175 IX4 = IX3+IDO
176 IF (NA .NE. 0) GO TO 110
177 CALL RADB5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
178 GO TO 111
179 110 CALL RADB5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
180 111 NA = 1-NA
181 GO TO 115
182 112 IF (NA .NE. 0) GO TO 113
183 CALL RADBG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
184 GO TO 114
185 113 CALL RADBG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
186 114 IF (IDO .EQ. 1) NA = 1-NA
187 115 L1 = L2
188 IW = IW+(IP-1)*IDO
189 116 CONTINUE
190 IF (NA .EQ. 0) RETURN
191 DO 117 I=1,N
192 C(I) = CH(I)
193 117 CONTINUE
194 RETURN
195 END
196 SUBROUTINE RFFTF1 (N,C,CH,WA,IFAC)
197 DIMENSION CH(1) ,C(1) ,WA(1) ,IFAC(*)
198 NF = IFAC(2)
199 NA = 1
200 L2 = N
201 IW = N
202 DO 111 K1=1,NF
203 KH = NF-K1
204 IP = IFAC(KH+3)
205 L1 = L2/IP
206 IDO = N/L2
207 IDL1 = IDO*L1
208 IW = IW-(IP-1)*IDO
209 NA = 1-NA
210 IF (IP .NE. 4) GO TO 102
211 IX2 = IW+IDO
212 IX3 = IX2+IDO
213 IF (NA .NE. 0) GO TO 101
214 CALL RADF4 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
215 GO TO 110
216 101 CALL RADF4 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
217 GO TO 110
218 102 IF (IP .NE. 2) GO TO 104
219 IF (NA .NE. 0) GO TO 103
220 CALL RADF2 (IDO,L1,C,CH,WA(IW))
221 GO TO 110
222 103 CALL RADF2 (IDO,L1,CH,C,WA(IW))
223 GO TO 110
224 104 IF (IP .NE. 3) GO TO 106
225 IX2 = IW+IDO
226 IF (NA .NE. 0) GO TO 105
227 CALL RADF3 (IDO,L1,C,CH,WA(IW),WA(IX2))
228 GO TO 110
229 105 CALL RADF3 (IDO,L1,CH,C,WA(IW),WA(IX2))
230 GO TO 110
231 106 IF (IP .NE. 5) GO TO 108
232 IX2 = IW+IDO
233 IX3 = IX2+IDO
234 IX4 = IX3+IDO
235 IF (NA .NE. 0) GO TO 107
236 CALL RADF5 (IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
237 GO TO 110
238 107 CALL RADF5 (IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
239 GO TO 110
240 108 IF (IDO .EQ. 1) NA = 1-NA
241 IF (NA .NE. 0) GO TO 109
242 CALL RADFG (IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
243 NA = 1
244 GO TO 110
245 109 CALL RADFG (IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
246 NA = 0
247 110 L2 = L1
248 111 CONTINUE
249 IF (NA .EQ. 1) RETURN
250 DO 112 I=1,N
251 C(I) = CH(I)
252 112 CONTINUE
253 RETURN
254 END
255 SUBROUTINE RFFTI1 (N,WA,IFAC)
256 DIMENSION WA(1) ,IFAC(*) ,NTRYH(4)
257 DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
258 NL = N
259 NF = 0
260 J = 0
261 101 J = J+1
262 IF (J-4) 102,102,103
263 102 NTRY = NTRYH(J)
264 GO TO 104
265 103 NTRY = NTRY+2
266 104 NQ = NL/NTRY
267 NR = NL-NTRY*NQ
268 IF (NR) 101,105,101
269 105 NF = NF+1
270 IFAC(NF+2) = NTRY
271 NL = NQ
272 IF (NTRY .NE. 2) GO TO 107
273 IF (NF .EQ. 1) GO TO 107
274 DO 106 I=2,NF
275 IB = NF-I+2
276 IFAC(IB+2) = IFAC(IB+1)
277 106 CONTINUE
278 IFAC(3) = 2
279 107 IF (NL .NE. 1) GO TO 104
280 IFAC(1) = N
281 IFAC(2) = NF
282 TPI = 6.28318530717959
283 ARGH = TPI/FLOAT(N)
284 IS = 0
285 NFM1 = NF-1
286 L1 = 1
287 IF (NFM1 .EQ. 0) RETURN
288
289 DO 110 K1=1,NFM1
290 IP = IFAC(K1+2)
291 LD = 0
292 L2 = L1*IP
293 IDO = N/L2
294 IPM = IP-1
295 DO 109 J=1,IPM
296 LD = LD+L1
297 I = IS
298 ARGLD = FLOAT(LD)*ARGH
299 FI = 0
300
301 DO 108 II=3,IDO,2
302 I = I+2
303 FI = FI+1
304 ARG = FI*ARGLD
305 WA(I-1) = COS(ARG)
306 WA(I) = SIN(ARG)
307 108 CONTINUE
308 IS = IS+IDO
309 109 CONTINUE
310 L1 = L2
311 110 CONTINUE
312 RETURN
313 END
314 SUBROUTINE RADB2 (IDO,L1,CC,CH,WA1)
315 DIMENSION CC(IDO,2,L1) ,CH(IDO,L1,2) ,
316 1 WA1(1)
317 DO 101 K=1,L1
318 CH(1,K,1) = CC(1,1,K)+CC(IDO,2,K)
319 CH(1,K,2) = CC(1,1,K)-CC(IDO,2,K)
320 101 CONTINUE
321 IF (IDO-2) 107,105,102
322 102 IDP2 = IDO+2
323
324 DO 104 K=1,L1
325 DO 103 I=3,IDO,2
326 IC = IDP2-I
327 CH(I-1,K,1) = CC(I-1,1,K)+CC(IC-1,2,K)
328 TR2 = CC(I-1,1,K)-CC(IC-1,2,K)
329 CH(I,K,1) = CC(I,1,K)-CC(IC,2,K)
330 TI2 = CC(I,1,K)+CC(IC,2,K)
331 CH(I-1,K,2) = WA1(I-2)*TR2-WA1(I-1)*TI2
332 CH(I,K,2) = WA1(I-2)*TI2+WA1(I-1)*TR2
333 103 CONTINUE
334 104 CONTINUE
335 IF (MOD(IDO,2) .EQ. 1) RETURN
336 105 DO 106 K=1,L1
337 CH(IDO,K,1) = CC(IDO,1,K)+CC(IDO,1,K)
338 CH(IDO,K,2) = -(CC(1,2,K)+CC(1,2,K))
339 106 CONTINUE
340 107 RETURN
341 END
342 SUBROUTINE RADB3 (IDO,L1,CC,CH,WA1,WA2)
343 DIMENSION CC(IDO,3,L1) ,CH(IDO,L1,3) ,
344 1 WA1(1) ,WA2(1)
345 DATA TAUR,TAUI /-.5,.866025403784439/
346 DO 101 K=1,L1
347 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
348 CR2 = CC(1,1,K)+TAUR*TR2
349 CH(1,K,1) = CC(1,1,K)+TR2
350 CI3 = TAUI*(CC(1,3,K)+CC(1,3,K))
351 CH(1,K,2) = CR2-CI3
352 CH(1,K,3) = CR2+CI3
353 101 CONTINUE
354 IF (IDO .EQ. 1) RETURN
355 IDP2 = IDO+2
356
357 DO 103 K=1,L1
358 DO 102 I=3,IDO,2
359 IC = IDP2-I
360 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
361 CR2 = CC(I-1,1,K)+TAUR*TR2
362 CH(I-1,K,1) = CC(I-1,1,K)+TR2
363 TI2 = CC(I,3,K)-CC(IC,2,K)
364 CI2 = CC(I,1,K)+TAUR*TI2
365 CH(I,K,1) = CC(I,1,K)+TI2
366 CR3 = TAUI*(CC(I-1,3,K)-CC(IC-1,2,K))
367 CI3 = TAUI*(CC(I,3,K)+CC(IC,2,K))
368 DR2 = CR2-CI3
369 DR3 = CR2+CI3
370 DI2 = CI2+CR3
371 DI3 = CI2-CR3
372 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
373 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
374 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
375 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
376 102 CONTINUE
377 103 CONTINUE
378 RETURN
379 END
380 SUBROUTINE RADB4 (IDO,L1,CC,CH,WA1,WA2,WA3)
381 DIMENSION CC(IDO,4,L1) ,CH(IDO,L1,4) ,
382 1 WA1(1) ,WA2(1) ,WA3(1)
383 DATA SQRT2 /1.414213562373095/
384 DO 101 K=1,L1
385 TR1 = CC(1,1,K)-CC(IDO,4,K)
386 TR2 = CC(1,1,K)+CC(IDO,4,K)
387 TR3 = CC(IDO,2,K)+CC(IDO,2,K)
388 TR4 = CC(1,3,K)+CC(1,3,K)
389 CH(1,K,1) = TR2+TR3
390 CH(1,K,2) = TR1-TR4
391 CH(1,K,3) = TR2-TR3
392 CH(1,K,4) = TR1+TR4
393 101 CONTINUE
394 IF (IDO-2) 107,105,102
395 102 IDP2 = IDO+2
396
397 DO 104 K=1,L1
398 DO 103 I=3,IDO,2
399 IC = IDP2-I
400 TI1 = CC(I,1,K)+CC(IC,4,K)
401 TI2 = CC(I,1,K)-CC(IC,4,K)
402 TI3 = CC(I,3,K)-CC(IC,2,K)
403 TR4 = CC(I,3,K)+CC(IC,2,K)
404 TR1 = CC(I-1,1,K)-CC(IC-1,4,K)
405 TR2 = CC(I-1,1,K)+CC(IC-1,4,K)
406 TI4 = CC(I-1,3,K)-CC(IC-1,2,K)
407 TR3 = CC(I-1,3,K)+CC(IC-1,2,K)
408 CH(I-1,K,1) = TR2+TR3
409 CR3 = TR2-TR3
410 CH(I,K,1) = TI2+TI3
411 CI3 = TI2-TI3
412 CR2 = TR1-TR4
413 CR4 = TR1+TR4
414 CI2 = TI1+TI4
415 CI4 = TI1-TI4
416 CH(I-1,K,2) = WA1(I-2)*CR2-WA1(I-1)*CI2
417 CH(I,K,2) = WA1(I-2)*CI2+WA1(I-1)*CR2
418 CH(I-1,K,3) = WA2(I-2)*CR3-WA2(I-1)*CI3
419 CH(I,K,3) = WA2(I-2)*CI3+WA2(I-1)*CR3
420 CH(I-1,K,4) = WA3(I-2)*CR4-WA3(I-1)*CI4
421 CH(I,K,4) = WA3(I-2)*CI4+WA3(I-1)*CR4
422 103 CONTINUE
423 104 CONTINUE
424 IF (MOD(IDO,2) .EQ. 1) RETURN
425 105 CONTINUE
426 DO 106 K=1,L1
427 TI1 = CC(1,2,K)+CC(1,4,K)
428 TI2 = CC(1,4,K)-CC(1,2,K)
429 TR1 = CC(IDO,1,K)-CC(IDO,3,K)
430 TR2 = CC(IDO,1,K)+CC(IDO,3,K)
431 CH(IDO,K,1) = TR2+TR2
432 CH(IDO,K,2) = SQRT2*(TR1-TI1)
433 CH(IDO,K,3) = TI2+TI2
434 CH(IDO,K,4) = -SQRT2*(TR1+TI1)
435 106 CONTINUE
436 107 RETURN
437 END
438 SUBROUTINE RADB5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
439 DIMENSION CC(IDO,5,L1) ,CH(IDO,L1,5) ,
440 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
441 DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
442 1-.809016994374947,.587785252292473/
443 DO 101 K=1,L1
444 TI5 = CC(1,3,K)+CC(1,3,K)
445 TI4 = CC(1,5,K)+CC(1,5,K)
446 TR2 = CC(IDO,2,K)+CC(IDO,2,K)
447 TR3 = CC(IDO,4,K)+CC(IDO,4,K)
448 CH(1,K,1) = CC(1,1,K)+TR2+TR3
449 CR2 = CC(1,1,K)+TR11*TR2+TR12*TR3
450 CR3 = CC(1,1,K)+TR12*TR2+TR11*TR3
451 CI5 = TI11*TI5+TI12*TI4
452 CI4 = TI12*TI5-TI11*TI4
453 CH(1,K,2) = CR2-CI5
454 CH(1,K,3) = CR3-CI4
455 CH(1,K,4) = CR3+CI4
456 CH(1,K,5) = CR2+CI5
457 101 CONTINUE
458 IF (IDO .EQ. 1) RETURN
459 IDP2 = IDO+2
460 DO 103 K=1,L1
461 DO 102 I=3,IDO,2
462 IC = IDP2-I
463 TI5 = CC(I,3,K)+CC(IC,2,K)
464 TI2 = CC(I,3,K)-CC(IC,2,K)
465 TI4 = CC(I,5,K)+CC(IC,4,K)
466 TI3 = CC(I,5,K)-CC(IC,4,K)
467 TR5 = CC(I-1,3,K)-CC(IC-1,2,K)
468 TR2 = CC(I-1,3,K)+CC(IC-1,2,K)
469 TR4 = CC(I-1,5,K)-CC(IC-1,4,K)
470 TR3 = CC(I-1,5,K)+CC(IC-1,4,K)
471 CH(I-1,K,1) = CC(I-1,1,K)+TR2+TR3
472 CH(I,K,1) = CC(I,1,K)+TI2+TI3
473 CR2 = CC(I-1,1,K)+TR11*TR2+TR12*TR3
474 CI2 = CC(I,1,K)+TR11*TI2+TR12*TI3
475 CR3 = CC(I-1,1,K)+TR12*TR2+TR11*TR3
476 CI3 = CC(I,1,K)+TR12*TI2+TR11*TI3
477 CR5 = TI11*TR5+TI12*TR4
478 CI5 = TI11*TI5+TI12*TI4
479 CR4 = TI12*TR5-TI11*TR4
480 CI4 = TI12*TI5-TI11*TI4
481 DR3 = CR3-CI4
482 DR4 = CR3+CI4
483 DI3 = CI3+CR4
484 DI4 = CI3-CR4
485 DR5 = CR2+CI5
486 DR2 = CR2-CI5
487 DI5 = CI2-CR5
488 DI2 = CI2+CR5
489 CH(I-1,K,2) = WA1(I-2)*DR2-WA1(I-1)*DI2
490 CH(I,K,2) = WA1(I-2)*DI2+WA1(I-1)*DR2
491 CH(I-1,K,3) = WA2(I-2)*DR3-WA2(I-1)*DI3
492 CH(I,K,3) = WA2(I-2)*DI3+WA2(I-1)*DR3
493 CH(I-1,K,4) = WA3(I-2)*DR4-WA3(I-1)*DI4
494 CH(I,K,4) = WA3(I-2)*DI4+WA3(I-1)*DR4
495 CH(I-1,K,5) = WA4(I-2)*DR5-WA4(I-1)*DI5
496 CH(I,K,5) = WA4(I-2)*DI5+WA4(I-1)*DR5
497 102 CONTINUE
498 103 CONTINUE
499 RETURN
500 END
501 SUBROUTINE RADBG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
502 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
503 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
504 2 CH2(IDL1,IP) ,WA(1)
505 DATA TPI/6.28318530717959/
506 ARG = TPI/FLOAT(IP)
507 DCP = COS(ARG)
508 DSP = SIN(ARG)
509 IDP2 = IDO+2
510 NBD = (IDO-1)/2
511 IPP2 = IP+2
512 IPPH = (IP+1)/2
513 IF (IDO .LT. L1) GO TO 103
514 DO 102 K=1,L1
515 DO 101 I=1,IDO
516 CH(I,K,1) = CC(I,1,K)
517 101 CONTINUE
518 102 CONTINUE
519 GO TO 106
520 103 DO 105 I=1,IDO
521 DO 104 K=1,L1
522 CH(I,K,1) = CC(I,1,K)
523 104 CONTINUE
524 105 CONTINUE
525
526 DO 108 J=2,IPPH
527 JC = IPP2-J
528 J2 = J+J
529 DO 107 K=1,L1
530 CH(1,K,J) = CC(IDO,J2-2,K)+CC(IDO,J2-2,K)
531 CH(1,K,JC) = CC(1,J2-1,K)+CC(1,J2-1,K)
532 107 CONTINUE
533 108 CONTINUE
534 IF (IDO .EQ. 1) GO TO 116
535 IF (NBD .LT. L1) GO TO 112
536
537 DO 111 J=2,IPPH
538 JC = IPP2-J
539 DO 110 K=1,L1
540 DO 109 I=3,IDO,2
541 IC = IDP2-I
542 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
543 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
544 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
545 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
546 109 CONTINUE
547 110 CONTINUE
548 111 CONTINUE
549 GO TO 116
550 112 DO 115 J=2,IPPH
551 JC = IPP2-J
552 DO 114 I=3,IDO,2
553 IC = IDP2-I
554 DO 113 K=1,L1
555 CH(I-1,K,J) = CC(I-1,2*J-1,K)+CC(IC-1,2*J-2,K)
556 CH(I-1,K,JC) = CC(I-1,2*J-1,K)-CC(IC-1,2*J-2,K)
557 CH(I,K,J) = CC(I,2*J-1,K)-CC(IC,2*J-2,K)
558 CH(I,K,JC) = CC(I,2*J-1,K)+CC(IC,2*J-2,K)
559 113 CONTINUE
560 114 CONTINUE
561 115 CONTINUE
562 116 AR1 = 1.
563 AI1 = 0.
564
565 DO 120 L=2,IPPH
566 LC = IPP2-L
567 AR1H = DCP*AR1-DSP*AI1
568 AI1 = DCP*AI1+DSP*AR1
569 AR1 = AR1H
570 DO 117 IK=1,IDL1
571 C2(IK,L) = CH2(IK,1)+AR1*CH2(IK,2)
572 C2(IK,LC) = AI1*CH2(IK,IP)
573 117 CONTINUE
574 DC2 = AR1
575 DS2 = AI1
576 AR2 = AR1
577 AI2 = AI1
578
579 DO 119 J=3,IPPH
580 JC = IPP2-J
581 AR2H = DC2*AR2-DS2*AI2
582 AI2 = DC2*AI2+DS2*AR2
583 AR2 = AR2H
584 DO 118 IK=1,IDL1
585 C2(IK,L) = C2(IK,L)+AR2*CH2(IK,J)
586 C2(IK,LC) = C2(IK,LC)+AI2*CH2(IK,JC)
587 118 CONTINUE
588 119 CONTINUE
589 120 CONTINUE
590
591 DO 122 J=2,IPPH
592 DO 121 IK=1,IDL1
593 CH2(IK,1) = CH2(IK,1)+CH2(IK,J)
594 121 CONTINUE
595 122 CONTINUE
596
597 DO 124 J=2,IPPH
598 JC = IPP2-J
599 DO 123 K=1,L1
600 CH(1,K,J) = C1(1,K,J)-C1(1,K,JC)
601 CH(1,K,JC) = C1(1,K,J)+C1(1,K,JC)
602 123 CONTINUE
603 124 CONTINUE
604 IF (IDO .EQ. 1) GO TO 132
605 IF (NBD .LT. L1) GO TO 128
606
607 DO 127 J=2,IPPH
608 JC = IPP2-J
609 DO 126 K=1,L1
610 DO 125 I=3,IDO,2
611 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
612 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
613 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
614 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
615 125 CONTINUE
616 126 CONTINUE
617 127 CONTINUE
618 GO TO 132
619 128 DO 131 J=2,IPPH
620 JC = IPP2-J
621 DO 130 I=3,IDO,2
622 DO 129 K=1,L1
623 CH(I-1,K,J) = C1(I-1,K,J)-C1(I,K,JC)
624 CH(I-1,K,JC) = C1(I-1,K,J)+C1(I,K,JC)
625 CH(I,K,J) = C1(I,K,J)+C1(I-1,K,JC)
626 CH(I,K,JC) = C1(I,K,J)-C1(I-1,K,JC)
627 129 CONTINUE
628 130 CONTINUE
629 131 CONTINUE
630 132 CONTINUE
631 IF (IDO .EQ. 1) RETURN
632 DO 133 IK=1,IDL1
633 C2(IK,1) = CH2(IK,1)
634 133 CONTINUE
635 DO 135 J=2,IP
636 DO 134 K=1,L1
637 C1(1,K,J) = CH(1,K,J)
638 134 CONTINUE
639 135 CONTINUE
640 IF (NBD .GT. L1) GO TO 139
641 IS = -IDO
642 DO 138 J=2,IP
643 IS = IS+IDO
644 IDIJ = IS
645 DO 137 I=3,IDO,2
646 IDIJ = IDIJ+2
647 DO 136 K=1,L1
648 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
649 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
650 136 CONTINUE
651 137 CONTINUE
652 138 CONTINUE
653 GO TO 143
654 139 IS = -IDO
655
656 DO 142 J=2,IP
657 IS = IS+IDO
658 DO 141 K=1,L1
659 IDIJ = IS
660 DO 140 I=3,IDO,2
661 IDIJ = IDIJ+2
662 C1(I-1,K,J) = WA(IDIJ-1)*CH(I-1,K,J)-WA(IDIJ)*CH(I,K,J)
663 C1(I,K,J) = WA(IDIJ-1)*CH(I,K,J)+WA(IDIJ)*CH(I-1,K,J)
664 140 CONTINUE
665 141 CONTINUE
666 142 CONTINUE
667 143 RETURN
668 END
669 SUBROUTINE RADF2 (IDO,L1,CC,CH,WA1)
670 DIMENSION CH(IDO,2,L1) ,CC(IDO,L1,2) ,
671 1 WA1(1)
672 DO 101 K=1,L1
673 CH(1,1,K) = CC(1,K,1)+CC(1,K,2)
674 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,2)
675 101 CONTINUE
676 IF (IDO-2) 107,105,102
677 102 IDP2 = IDO+2
678 DO 104 K=1,L1
679 DO 103 I=3,IDO,2
680 IC = IDP2-I
681 TR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
682 TI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
683 CH(I,1,K) = CC(I,K,1)+TI2
684 CH(IC,2,K) = TI2-CC(I,K,1)
685 CH(I-1,1,K) = CC(I-1,K,1)+TR2
686 CH(IC-1,2,K) = CC(I-1,K,1)-TR2
687 103 CONTINUE
688 104 CONTINUE
689 IF (MOD(IDO,2) .EQ. 1) RETURN
690 105 DO 106 K=1,L1
691 CH(1,2,K) = -CC(IDO,K,2)
692 CH(IDO,1,K) = CC(IDO,K,1)
693 106 CONTINUE
694 107 RETURN
695 END
696 SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)
697 DIMENSION CH(IDO,3,L1) ,CC(IDO,L1,3) ,
698 1 WA1(1) ,WA2(1)
699 DATA TAUR,TAUI /-.5,.866025403784439/
700 DO 101 K=1,L1
701 CR2 = CC(1,K,2)+CC(1,K,3)
702 CH(1,1,K) = CC(1,K,1)+CR2
703 CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
704 CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
705 101 CONTINUE
706 IF (IDO .EQ. 1) RETURN
707 IDP2 = IDO+2
708 DO 103 K=1,L1
709 DO 102 I=3,IDO,2
710 IC = IDP2-I
711 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
712 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
713 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
714 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
715 CR2 = DR2+DR3
716 CI2 = DI2+DI3
717 CH(I-1,1,K) = CC(I-1,K,1)+CR2
718 CH(I,1,K) = CC(I,K,1)+CI2
719 TR2 = CC(I-1,K,1)+TAUR*CR2
720 TI2 = CC(I,K,1)+TAUR*CI2
721 TR3 = TAUI*(DI2-DI3)
722 TI3 = TAUI*(DR3-DR2)
723 CH(I-1,3,K) = TR2+TR3
724 CH(IC-1,2,K) = TR2-TR3
725 CH(I,3,K) = TI2+TI3
726 CH(IC,2,K) = TI3-TI2
727 102 CONTINUE
728 103 CONTINUE
729 RETURN
730 END
731 SUBROUTINE RADF4 (IDO,L1,CC,CH,WA1,WA2,WA3)
732 DIMENSION CC(IDO,L1,4) ,CH(IDO,4,L1) ,
733 1 WA1(1) ,WA2(1) ,WA3(1)
734 DATA HSQT2 /.7071067811865475/
735 DO 101 K=1,L1
736 TR1 = CC(1,K,2)+CC(1,K,4)
737 TR2 = CC(1,K,1)+CC(1,K,3)
738 CH(1,1,K) = TR1+TR2
739 CH(IDO,4,K) = TR2-TR1
740 CH(IDO,2,K) = CC(1,K,1)-CC(1,K,3)
741 CH(1,3,K) = CC(1,K,4)-CC(1,K,2)
742 101 CONTINUE
743 IF (IDO-2) 107,105,102
744 102 IDP2 = IDO+2
745
746 DO 104 K=1,L1
747 DO 103 I=3,IDO,2
748 IC = IDP2-I
749 CR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
750 CI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
751 CR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
752 CI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
753 CR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
754 CI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
755 TR1 = CR2+CR4
756 TR4 = CR4-CR2
757 TI1 = CI2+CI4
758 TI4 = CI2-CI4
759 TI2 = CC(I,K,1)+CI3
760 TI3 = CC(I,K,1)-CI3
761 TR2 = CC(I-1,K,1)+CR3
762 TR3 = CC(I-1,K,1)-CR3
763 CH(I-1,1,K) = TR1+TR2
764 CH(IC-1,4,K) = TR2-TR1
765 CH(I,1,K) = TI1+TI2
766 CH(IC,4,K) = TI1-TI2
767 CH(I-1,3,K) = TI4+TR3
768 CH(IC-1,2,K) = TR3-TI4
769 CH(I,3,K) = TR4+TI3
770 CH(IC,2,K) = TR4-TI3
771 103 CONTINUE
772 104 CONTINUE
773 IF (MOD(IDO,2) .EQ. 1) RETURN
774 105 CONTINUE
775 DO 106 K=1,L1
776 TI1 = -HSQT2*(CC(IDO,K,2)+CC(IDO,K,4))
777 TR1 = HSQT2*(CC(IDO,K,2)-CC(IDO,K,4))
778 CH(IDO,1,K) = TR1+CC(IDO,K,1)
779 CH(IDO,3,K) = CC(IDO,K,1)-TR1
780 CH(1,2,K) = TI1-CC(IDO,K,3)
781 CH(1,4,K) = TI1+CC(IDO,K,3)
782 106 CONTINUE
783 107 RETURN
784 END
785 SUBROUTINE RADF5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
786 DIMENSION CC(IDO,L1,5) ,CH(IDO,5,L1) ,
787 1 WA1(1) ,WA2(1) ,WA3(1) ,WA4(1)
788 DATA TR11,TI11,TR12,TI12 /.309016994374947,.951056516295154,
789 1-.809016994374947,.587785252292473/
790 DO 101 K=1,L1
791 CR2 = CC(1,K,5)+CC(1,K,2)
792 CI5 = CC(1,K,5)-CC(1,K,2)
793 CR3 = CC(1,K,4)+CC(1,K,3)
794 CI4 = CC(1,K,4)-CC(1,K,3)
795 CH(1,1,K) = CC(1,K,1)+CR2+CR3
796 CH(IDO,2,K) = CC(1,K,1)+TR11*CR2+TR12*CR3
797 CH(1,3,K) = TI11*CI5+TI12*CI4
798 CH(IDO,4,K) = CC(1,K,1)+TR12*CR2+TR11*CR3
799 CH(1,5,K) = TI12*CI5-TI11*CI4
800 101 CONTINUE
801 IF (IDO .EQ. 1) RETURN
802 IDP2 = IDO+2
803 DO 103 K=1,L1
804 DO 102 I=3,IDO,2
805 IC = IDP2-I
806 DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
807 DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
808 DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
809 DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
810 DR4 = WA3(I-2)*CC(I-1,K,4)+WA3(I-1)*CC(I,K,4)
811 DI4 = WA3(I-2)*CC(I,K,4)-WA3(I-1)*CC(I-1,K,4)
812 DR5 = WA4(I-2)*CC(I-1,K,5)+WA4(I-1)*CC(I,K,5)
813 DI5 = WA4(I-2)*CC(I,K,5)-WA4(I-1)*CC(I-1,K,5)
814 CR2 = DR2+DR5
815 CI5 = DR5-DR2
816 CR5 = DI2-DI5
817 CI2 = DI2+DI5
818 CR3 = DR3+DR4
819 CI4 = DR4-DR3
820 CR4 = DI3-DI4
821 CI3 = DI3+DI4
822 CH(I-1,1,K) = CC(I-1,K,1)+CR2+CR3
823 CH(I,1,K) = CC(I,K,1)+CI2+CI3
824 TR2 = CC(I-1,K,1)+TR11*CR2+TR12*CR3
825 TI2 = CC(I,K,1)+TR11*CI2+TR12*CI3
826 TR3 = CC(I-1,K,1)+TR12*CR2+TR11*CR3
827 TI3 = CC(I,K,1)+TR12*CI2+TR11*CI3
828 TR5 = TI11*CR5+TI12*CR4
829 TI5 = TI11*CI5+TI12*CI4
830 TR4 = TI12*CR5-TI11*CR4
831 TI4 = TI12*CI5-TI11*CI4
832 CH(I-1,3,K) = TR2+TR5
833 CH(IC-1,2,K) = TR2-TR5
834 CH(I,3,K) = TI2+TI5
835 CH(IC,2,K) = TI5-TI2
836 CH(I-1,5,K) = TR3+TR4
837 CH(IC-1,4,K) = TR3-TR4
838 CH(I,5,K) = TI3+TI4
839 CH(IC,4,K) = TI4-TI3
840 102 CONTINUE
841 103 CONTINUE
842 RETURN
843 END
844 SUBROUTINE RADFG (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
845 DIMENSION CH(IDO,L1,IP) ,CC(IDO,IP,L1) ,
846 1 C1(IDO,L1,IP) ,C2(IDL1,IP),
847 2 CH2(IDL1,IP) ,WA(1)
848 DATA TPI/6.28318530717959/
849 ARG = TPI/FLOAT(IP)
850 DCP = COS(ARG)
851 DSP = SIN(ARG)
852 IPPH = (IP+1)/2
853 IPP2 = IP+2
854 IDP2 = IDO+2
855 NBD = (IDO-1)/2
856 IF (IDO .EQ. 1) GO TO 119
857 DO 101 IK=1,IDL1
858 CH2(IK,1) = C2(IK,1)
859 101 CONTINUE
860 DO 103 J=2,IP
861 DO 102 K=1,L1
862 CH(1,K,J) = C1(1,K,J)
863 102 CONTINUE
864 103 CONTINUE
865 IF (NBD .GT. L1) GO TO 107
866 IS = -IDO
867 DO 106 J=2,IP
868 IS = IS+IDO
869 IDIJ = IS
870 DO 105 I=3,IDO,2
871 IDIJ = IDIJ+2
872 DO 104 K=1,L1
873 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
874 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
875 104 CONTINUE
876 105 CONTINUE
877 106 CONTINUE
878 GO TO 111
879 107 IS = -IDO
880 DO 110 J=2,IP
881 IS = IS+IDO
882 DO 109 K=1,L1
883 IDIJ = IS
884 DO 108 I=3,IDO,2
885 IDIJ = IDIJ+2
886 CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J)+WA(IDIJ)*C1(I,K,J)
887 CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J)-WA(IDIJ)*C1(I-1,K,J)
888 108 CONTINUE
889 109 CONTINUE
890 110 CONTINUE
891 111 IF (NBD .LT. L1) GO TO 115
892 DO 114 J=2,IPPH
893 JC = IPP2-J
894 DO 113 K=1,L1
895 DO 112 I=3,IDO,2
896 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
897 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
898 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
899 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
900 112 CONTINUE
901 113 CONTINUE
902 114 CONTINUE
903 GO TO 121
904 115 DO 118 J=2,IPPH
905 JC = IPP2-J
906 DO 117 I=3,IDO,2
907 DO 116 K=1,L1
908 C1(I-1,K,J) = CH(I-1,K,J)+CH(I-1,K,JC)
909 C1(I-1,K,JC) = CH(I,K,J)-CH(I,K,JC)
910 C1(I,K,J) = CH(I,K,J)+CH(I,K,JC)
911 C1(I,K,JC) = CH(I-1,K,JC)-CH(I-1,K,J)
912 116 CONTINUE
913 117 CONTINUE
914 118 CONTINUE
915 GO TO 121
916 119 DO 120 IK=1,IDL1
917 C2(IK,1) = CH2(IK,1)
918 120 CONTINUE
919 121 DO 123 J=2,IPPH
920 JC = IPP2-J
921 DO 122 K=1,L1
922 C1(1,K,J) = CH(1,K,J)+CH(1,K,JC)
923 C1(1,K,JC) = CH(1,K,JC)-CH(1,K,J)
924 122 CONTINUE
925 123 CONTINUE
926
927 = 1.
928 AI1 = 0.
929 DO 127 L=2,IPPH
930 LC = IPP2-L
931 AR1H = DCP*AR1-DSP*AI1
932 AI1 = DCP*AI1+DSP*AR1
933 AR1 = AR1H
934 DO 124 IK=1,IDL1
935 CH2(IK,L) = C2(IK,1)+AR1*C2(IK,2)
936 CH2(IK,LC) = AI1*C2(IK,IP)
937 124 CONTINUE
938 DC2 = AR1
939 DS2 = AI1
940 AR2 = AR1
941 AI2 = AI1
942 DO 126 J=3,IPPH
943 JC = IPP2-J
944 AR2H = DC2*AR2-DS2*AI2
945 AI2 = DC2*AI2+DS2*AR2
946 AR2 = AR2H
947 DO 125 IK=1,IDL1
948 CH2(IK,L) = CH2(IK,L)+AR2*C2(IK,J)
949 CH2(IK,LC) = CH2(IK,LC)+AI2*C2(IK,JC)
950 125 CONTINUE
951 126 CONTINUE
952 127 CONTINUE
953 DO 129 J=2,IPPH
954 DO 128 IK=1,IDL1
955 CH2(IK,1) = CH2(IK,1)+C2(IK,J)
956 128 CONTINUE
957 129 CONTINUE
958
959 IF (IDO .LT. L1) GO TO 132
960 DO 131 K=1,L1
961 DO 130 I=1,IDO
962 CC(I,1,K) = CH(I,K,1)
963 130 CONTINUE
964 131 CONTINUE
965 GO TO 135
966 132 DO 134 I=1,IDO
967 DO 133 K=1,L1
968 CC(I,1,K) = CH(I,K,1)
969 133 CONTINUE
970 134 CONTINUE
971 135 DO 137 J=2,IPPH
972 JC = IPP2-J
973 J2 = J+J
974 DO 136 K=1,L1
975 CC(IDO,J2-2,K) = CH(1,K,J)
976 CC(1,J2-1,K) = CH(1,K,JC)
977 136 CONTINUE
978 137 CONTINUE
979 IF (IDO .EQ. 1) RETURN
980 IF (NBD .LT. L1) GO TO 141
981 DO 140 J=2,IPPH
982 JC = IPP2-J
983 J2 = J+J
984 DO 139 K=1,L1
985 DO 138 I=3,IDO,2
986 IC = IDP2-I
987 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
988 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
989 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
990 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
991 138 CONTINUE
992 139 CONTINUE
993 140 CONTINUE
994 RETURN
995 141 DO 144 J=2,IPPH
996 JC = IPP2-J
997 J2 = J+J
998 DO 143 I=3,IDO,2
999 IC = IDP2-I
1000 DO 142 K=1,L1
1001 CC(I-1,J2-1,K) = CH(I-1,K,J)+CH(I-1,K,JC)
1002 CC(IC-1,J2-2,K) = CH(I-1,K,J)-CH(I-1,K,JC)
1003 CC(I,J2-1,K) = CH(I,K,J)+CH(I,K,JC)
1004 CC(IC,J2-2,K) = CH(I,K,JC)-CH(I,K,J)
1005 142 CONTINUE
1006 143 CONTINUE
1007 144 CONTINUE
1008 RETURN
1009 END
1010