1*
2* $Id$
3*
4CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
5C
6C      File:           sint.f
7C
8C      Library:        FFTPACK 4.1
9C
10C      Author:         Paul N. Swarztrauber
11C                      National Center for Atmospheric Research
12C                      PO 3000, Boulder, Colorado
13C
14C      Date:           Wed Mar 29 18:31:13 MST 1995
15C
16C      Description:    Forward and backward, 1D real sin FFT
17C
18      SUBROUTINE SINT(N,X,WSAVE)
19      DOUBLE PRECISION X(*),WSAVE(*)
20
21      NP1 = N + 1
22      IW1 = N/2 + 1
23      IW2 = IW1 + NP1
24      IW3 = IW2 + NP1
25      CALL SINT1(N,X,WSAVE,WSAVE(IW1),WSAVE(IW2),WSAVE(IW3))
26      RETURN
27      END
28CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
29C
30C      File:           sinti.f
31C
32C      Library:        FFTPACK 4.1
33C
34C      Author:         Paul N. Swarztrauber
35C                      National Center for Atmospheric Research
36C                      PO 3000, Boulder, Colorado
37C
38C      Date:           Wed Mar 29 18:31:13 MST 1995
39C
40C      Description:    Initialization routine for SINT
41C
42      SUBROUTINE SINTI(N,WSAVE)
43      DOUBLE PRECISION PI
44      DOUBLE PRECISION PIMACH
45      DOUBLE PRECISION DUM
46      DOUBLE PRECISION DT
47      DOUBLE PRECISION WSAVE(*)
48
49      PI = PIMACH(DUM)
50      IF (N.LE.1) RETURN
51      NS2 = N/2
52      NP1 = N + 1
53      DT = PI/DBLE(NP1)
54      DO 101 K = 1,NS2
55          WSAVE(K) = 2.D0*SIN(K*DT)
56  101 CONTINUE
57      CALL RFFTI(NP1,WSAVE(NS2+1))
58      RETURN
59      END
60CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
61C
62C      File:           sint1.f
63C
64C      Library:        FFTPACK 4.1
65C
66C      Author:         Paul N. Swarztrauber
67C                      National Center for Atmospheric Research
68C                      PO 3000, Boulder, Colorado
69C
70C      Date:           Wed Mar 29 18:31:13 MST 1995
71C
72C      Description:    Lower-level auxiliary routine
73C
74      SUBROUTINE SINT1(N,WAR,WAS,XH,X,IFAC)
75      DOUBLE PRECISION SQRT3
76      DOUBLE PRECISION XHOLD
77      DOUBLE PRECISION T1
78      DOUBLE PRECISION T2
79      DOUBLE PRECISION WAR(*),WAS(*),X(*),XH(*)
80      integer*4 n4
81      INTEGER IFAC(*)
82      DATA SQRT3/1.73205080756888D0/
83C
84C FFTPACK 5.0 auxiliary routine
85C
86      DO 100 I = 1,N
87          XH(I) = WAR(I)
88          WAR(I) = X(I)
89  100 CONTINUE
90#if 1
91      n4=n-2
92      IF (n4) 101,102,103
93#else
94      IF (N-2) 101,102,103
95#endif
96  101 XH(1) = XH(1) + XH(1)
97      GO TO 106
98  102 XHOLD = SQRT3* (XH(1)+XH(2))
99      XH(2) = SQRT3* (XH(1)-XH(2))
100      XH(1) = XHOLD
101      GO TO 106
102  103 NP1 = N + 1
103      NS2 = N/2
104      X(1) = 0.D0
105      DO 104 K = 1,NS2
106          KC = NP1 - K
107          T1 = XH(K) - XH(KC)
108          T2 = WAS(K)* (XH(K)+XH(KC))
109          X(K+1) = T1 + T2
110          X(KC+1) = T2 - T1
111  104 CONTINUE
112      MODN = MOD(N,2)
113      IF (MODN.NE.0) X(NS2+2) = 4.D0*XH(NS2+1)
114      CALL RFFTF1(NP1,X,XH,WAR,IFAC)
115      XH(1) = .5D0*X(1)
116      DO 105 I = 3,N,2
117          XH(I-1) = -X(I)
118          XH(I) = XH(I-2) + X(I-1)
119  105 CONTINUE
120      IF (MODN.NE.0) GO TO 106
121      XH(N) = -X(N+1)
122  106 DO 107 I = 1,N
123          X(I) = WAR(I)
124          WAR(I) = XH(I)
125  107 CONTINUE
126      RETURN
127      END
128CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
129C
130C      File:           rffti.f
131C
132C      Library:        FFTPACK 4.1
133C
134C      Author:         Paul N. Swarztrauber
135C                      National Center for Atmospheric Research
136C                      PO 3000, Boulder, Colorado
137C
138C      Date:           Wed Mar 29 18:31:13 MST 1995
139C
140C      Description:    Initialization routine for RFFTB, RFFTF
141C
142      SUBROUTINE RFFTI(N,WSAVE)
143      DOUBLE PRECISION WSAVE(*)
144
145      IF (N.EQ.1) RETURN
146      CALL RFFTI1(N,WSAVE(N+1),WSAVE(2*N+1))
147      RETURN
148      END
149CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
150C
151C      File:           rfftf1.f
152C
153C      Library:        FFTPACK 4.1
154C
155C      Author:         Paul N. Swarztrauber
156C                      National Center for Atmospheric Research
157C                      PO 3000, Boulder, Colorado
158C
159C      Date:           Wed Mar 29 18:31:13 MST 1995
160C
161C      Description:    Lower-level auxiliary routine
162C
163      SUBROUTINE RFFTF1(N,C,CH,WA,IFAC)
164      DOUBLE PRECISION CH(*),C(*),WA(*)
165      INTEGER IFAC(*)
166C
167C FFTPACK 5.0 auxiliary routine
168C
169      NF = IFAC(2)
170      NA = 1
171      L2 = N
172      IW = N
173      DO 111 K1 = 1,NF
174          KH = NF - K1
175          IP = IFAC(KH+3)
176          L1 = L2/IP
177          IDO = N/L2
178          IDL1 = IDO*L1
179          IW = IW - (IP-1)*IDO
180          NA = 1 - NA
181          IF (IP.NE.4) GO TO 102
182          IX2 = IW + IDO
183          IX3 = IX2 + IDO
184          IF (NA.NE.0) GO TO 101
185          CALL RADF4(IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3))
186          GO TO 110
187  101     CALL RADF4(IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3))
188          GO TO 110
189  102     IF (IP.NE.2) GO TO 104
190          IF (NA.NE.0) GO TO 103
191          CALL RADF2(IDO,L1,C,CH,WA(IW))
192          GO TO 110
193  103     CALL RADF2(IDO,L1,CH,C,WA(IW))
194          GO TO 110
195  104     IF (IP.NE.3) GO TO 106
196          IX2 = IW + IDO
197          IF (NA.NE.0) GO TO 105
198          CALL RADF3(IDO,L1,C,CH,WA(IW),WA(IX2))
199          GO TO 110
200  105     CALL RADF3(IDO,L1,CH,C,WA(IW),WA(IX2))
201          GO TO 110
202  106     IF (IP.NE.5) GO TO 108
203          IX2 = IW + IDO
204          IX3 = IX2 + IDO
205          IX4 = IX3 + IDO
206          IF (NA.NE.0) GO TO 107
207          CALL RADF5(IDO,L1,C,CH,WA(IW),WA(IX2),WA(IX3),WA(IX4))
208          GO TO 110
209  107     CALL RADF5(IDO,L1,CH,C,WA(IW),WA(IX2),WA(IX3),WA(IX4))
210          GO TO 110
211  108     IF (IDO.EQ.1) NA = 1 - NA
212          IF (NA.NE.0) GO TO 109
213          CALL RADFG(IDO,IP,L1,IDL1,C,C,C,CH,CH,WA(IW))
214          NA = 1
215          GO TO 110
216  109     CALL RADFG(IDO,IP,L1,IDL1,CH,CH,CH,C,C,WA(IW))
217          NA = 0
218  110     L2 = L1
219  111 CONTINUE
220      IF (NA.EQ.1) RETURN
221      DO 112 I = 1,N
222          C(I) = CH(I)
223  112 CONTINUE
224      RETURN
225      END
226CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
227C
228C      File:           rffti1.f
229C
230C      Library:        FFTPACK 4.1
231C
232C      Author:         Paul N. Swarztrauber
233C                      National Center for Atmospheric Research
234C                      PO 3000, Boulder, Colorado
235C
236C      Date:           Wed Mar 29 18:31:13 MST 1995
237C
238C      Description:    Lower-level auxiliary routine
239C
240      SUBROUTINE RFFTI1(N,WA,IFAC)
241      DOUBLE PRECISION TPI
242      DOUBLE PRECISION PIMACH
243      DOUBLE PRECISION DUM
244      DOUBLE PRECISION ARGH
245      DOUBLE PRECISION ARGLD
246      DOUBLE PRECISION FI
247      DOUBLE PRECISION ARG
248      DOUBLE PRECISION WA(*)
249      INTEGER IFAC(*),NTRYH(4)
250      integer*4 n4
251      DATA NTRYH(1),NTRYH(2),NTRYH(3),NTRYH(4)/4,2,3,5/
252C
253C FFTPACK 5.0 auxiliary routine
254C
255      NL = N
256      NF = 0
257      NTRY = 0
258      J = 0
259  101 J = J + 1
260#if 1
261      n4=j-4
262      IF (n4) 102,102,103
263#else
264      IF (J-4) 102,102,103
265#endif
266  102 NTRY = NTRYH(J)
267      GO TO 104
268  103 NTRY = NTRY + 2
269  104 NQ = NL/NTRY
270      NR = NL - NTRY*NQ
271#if 1
272      n4=NR
273      IF (n4) 101,105,101
274#else
275      IF (NR) 101,105,101
276#endif
277  105 NF = NF + 1
278      IFAC(NF+2) = NTRY
279      NL = NQ
280      IF (NTRY.NE.2) GO TO 107
281      IF (NF.EQ.1) GO TO 107
282      DO 106 I = 2,NF
283          IB = NF - I + 2
284          IFAC(IB+2) = IFAC(IB+1)
285  106 CONTINUE
286      IFAC(3) = 2
287  107 IF (NL.NE.1) GO TO 104
288      IFAC(1) = N
289      IFAC(2) = NF
290      TPI = 2.0D0*PIMACH(DUM)
291      ARGH = TPI/DBLE(N)
292      IS = 0
293      NFM1 = NF - 1
294      L1 = 1
295      IF (NFM1.EQ.0) RETURN
296      DO 110 K1 = 1,NFM1
297          IP = IFAC(K1+2)
298          LD = 0
299          L2 = L1*IP
300          IDO = N/L2
301          IPM = IP - 1
302          DO 109 J = 1,IPM
303              LD = LD + L1
304              I = IS
305              ARGLD = DBLE(LD)*ARGH
306              FI = 0.D0
307              DO 108 II = 3,IDO,2
308                  I = I + 2
309                  FI = FI + 1.D0
310                  ARG = FI*ARGLD
311                  WA(I-1) = COS(ARG)
312                  WA(I) = SIN(ARG)
313  108         CONTINUE
314              IS = IS + IDO
315  109     CONTINUE
316          L1 = L2
317  110 CONTINUE
318      RETURN
319      END
320CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
321C
322C      File:           radf2.f
323C
324C      Library:        FFTPACK 4.1
325C
326C      Author:         Paul N. Swarztrauber
327C                      National Center for Atmospheric Research
328C                      PO 3000, Boulder, Colorado
329C
330C      Date:           Wed Mar 29 18:31:13 MST 1995
331C
332C      Description:    Lower-level auxiliary routine
333C
334      SUBROUTINE RADF2(IDO,L1,CC,CH,WA1)
335      DOUBLE PRECISION TR2
336      DOUBLE PRECISION TI2
337      DOUBLE PRECISION CH(IDO,2,L1),CC(IDO,L1,2),WA1(*)
338      integer*4 n4
339C
340C FFTPACK 5.0 auxiliary routine
341C
342      DO 101 K = 1,L1
343          CH(1,1,K) = CC(1,K,1) + CC(1,K,2)
344          CH(IDO,2,K) = CC(1,K,1) - CC(1,K,2)
345  101 CONTINUE
346#if 1
347      n4=ido-2
348      IF (n4) 107,105,102
349#else
350      IF (IDO-2) 107,105,102
351#endif
352  102 IDP2 = IDO + 2
353      DO 104 K = 1,L1
354          DO 103 I = 3,IDO,2
355              IC = IDP2 - I
356              TR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2)
357              TI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2)
358              CH(I,1,K) = CC(I,K,1) + TI2
359              CH(IC,2,K) = TI2 - CC(I,K,1)
360              CH(I-1,1,K) = CC(I-1,K,1) + TR2
361              CH(IC-1,2,K) = CC(I-1,K,1) - TR2
362  103     CONTINUE
363  104 CONTINUE
364      IF (MOD(IDO,2).EQ.1) RETURN
365  105 DO 106 K = 1,L1
366          CH(1,2,K) = -CC(IDO,K,2)
367          CH(IDO,1,K) = CC(IDO,K,1)
368  106 CONTINUE
369  107 RETURN
370      END
371CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
372C
373C      File:           radf3.f
374C
375C      Library:        FFTPACK 4.1
376C
377C      Author:         Paul N. Swarztrauber
378C                      National Center for Atmospheric Research
379C                      PO 3000, Boulder, Colorado
380C
381C      Date:           Wed Mar 29 18:31:13 MST 1995
382C
383C      Description:    Lower-level auxiliary routine
384C
385      SUBROUTINE RADF3(IDO,L1,CC,CH,WA1,WA2)
386      DOUBLE PRECISION TAUR
387      DOUBLE PRECISION TAUI
388      DOUBLE PRECISION CR2
389      DOUBLE PRECISION DR2
390      DOUBLE PRECISION DI2
391      DOUBLE PRECISION DR3
392      DOUBLE PRECISION DI3
393      DOUBLE PRECISION CI2
394      DOUBLE PRECISION TR2
395      DOUBLE PRECISION TI2
396      DOUBLE PRECISION TR3
397      DOUBLE PRECISION TI3
398      DOUBLE PRECISION CH(IDO,3,L1),CC(IDO,L1,3),WA1(*),WA2(*)
399      DATA TAUR,TAUI/-.5D0,.866025403784439D0/
400C
401C FFTPACK 5.0 auxiliary routine
402C
403      DO 101 K = 1,L1
404          CR2 = CC(1,K,2) + CC(1,K,3)
405          CH(1,1,K) = CC(1,K,1) + CR2
406          CH(1,3,K) = TAUI* (CC(1,K,3)-CC(1,K,2))
407          CH(IDO,2,K) = CC(1,K,1) + TAUR*CR2
408  101 CONTINUE
409      IF (IDO.EQ.1) RETURN
410      IDP2 = IDO + 2
411      DO 103 K = 1,L1
412          DO 102 I = 3,IDO,2
413              IC = IDP2 - I
414              DR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2)
415              DI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2)
416              DR3 = WA2(I-2)*CC(I-1,K,3) + WA2(I-1)*CC(I,K,3)
417              DI3 = WA2(I-2)*CC(I,K,3) - WA2(I-1)*CC(I-1,K,3)
418              CR2 = DR2 + DR3
419              CI2 = DI2 + DI3
420              CH(I-1,1,K) = CC(I-1,K,1) + CR2
421              CH(I,1,K) = CC(I,K,1) + CI2
422              TR2 = CC(I-1,K,1) + TAUR*CR2
423              TI2 = CC(I,K,1) + TAUR*CI2
424              TR3 = TAUI* (DI2-DI3)
425              TI3 = TAUI* (DR3-DR2)
426              CH(I-1,3,K) = TR2 + TR3
427              CH(IC-1,2,K) = TR2 - TR3
428              CH(I,3,K) = TI2 + TI3
429              CH(IC,2,K) = TI3 - TI2
430  102     CONTINUE
431  103 CONTINUE
432      RETURN
433      END
434CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
435C
436C      File:           radf4.f
437C
438C      Library:        FFTPACK 4.1
439C
440C      Author:         Paul N. Swarztrauber
441C                      National Center for Atmospheric Research
442C                      PO 3000, Boulder, Colorado
443C
444C      Date:           Wed Mar 29 18:31:13 MST 1995
445C
446C      Description:    Lower-level auxiliary routine
447C
448      SUBROUTINE RADF4(IDO,L1,CC,CH,WA1,WA2,WA3)
449      DOUBLE PRECISION HSQT2
450      DOUBLE PRECISION TR1
451      DOUBLE PRECISION TR2
452      DOUBLE PRECISION CR2
453      DOUBLE PRECISION CI2
454      DOUBLE PRECISION CR3
455      DOUBLE PRECISION CI3
456      DOUBLE PRECISION CR4
457      DOUBLE PRECISION CI4
458      DOUBLE PRECISION TR4
459      DOUBLE PRECISION TI1
460      DOUBLE PRECISION TI4
461      DOUBLE PRECISION TI2
462      DOUBLE PRECISION TI3
463      DOUBLE PRECISION TR3
464      DOUBLE PRECISION CC(IDO,L1,4),CH(IDO,4,L1),WA1(*),WA2(*),WA3(*)
465      integer*4 n4
466      DATA HSQT2/.7071067811865475D0/
467C
468C FFTPACK 5.0 auxiliary routine
469C
470      DO 101 K = 1,L1
471          TR1 = CC(1,K,2) + CC(1,K,4)
472          TR2 = CC(1,K,1) + CC(1,K,3)
473          CH(1,1,K) = TR1 + TR2
474          CH(IDO,4,K) = TR2 - TR1
475          CH(IDO,2,K) = CC(1,K,1) - CC(1,K,3)
476          CH(1,3,K) = CC(1,K,4) - CC(1,K,2)
477  101 CONTINUE
478#if 1
479      n4=ido-2
480      IF (n4) 107,105,102
481#else
482      IF (IDO-2) 107,105,102
483#endif
484  102 IDP2 = IDO + 2
485      DO 104 K = 1,L1
486          DO 103 I = 3,IDO,2
487              IC = IDP2 - I
488              CR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2)
489              CI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2)
490              CR3 = WA2(I-2)*CC(I-1,K,3) + WA2(I-1)*CC(I,K,3)
491              CI3 = WA2(I-2)*CC(I,K,3) - WA2(I-1)*CC(I-1,K,3)
492              CR4 = WA3(I-2)*CC(I-1,K,4) + WA3(I-1)*CC(I,K,4)
493              CI4 = WA3(I-2)*CC(I,K,4) - WA3(I-1)*CC(I-1,K,4)
494              TR1 = CR2 + CR4
495              TR4 = CR4 - CR2
496              TI1 = CI2 + CI4
497              TI4 = CI2 - CI4
498              TI2 = CC(I,K,1) + CI3
499              TI3 = CC(I,K,1) - CI3
500              TR2 = CC(I-1,K,1) + CR3
501              TR3 = CC(I-1,K,1) - CR3
502              CH(I-1,1,K) = TR1 + TR2
503              CH(IC-1,4,K) = TR2 - TR1
504              CH(I,1,K) = TI1 + TI2
505              CH(IC,4,K) = TI1 - TI2
506              CH(I-1,3,K) = TI4 + TR3
507              CH(IC-1,2,K) = TR3 - TI4
508              CH(I,3,K) = TR4 + TI3
509              CH(IC,2,K) = TR4 - TI3
510  103     CONTINUE
511  104 CONTINUE
512      IF (MOD(IDO,2).EQ.1) RETURN
513  105 CONTINUE
514      DO 106 K = 1,L1
515          TI1 = -HSQT2* (CC(IDO,K,2)+CC(IDO,K,4))
516          TR1 = HSQT2* (CC(IDO,K,2)-CC(IDO,K,4))
517          CH(IDO,1,K) = TR1 + CC(IDO,K,1)
518          CH(IDO,3,K) = CC(IDO,K,1) - TR1
519          CH(1,2,K) = TI1 - CC(IDO,K,3)
520          CH(1,4,K) = TI1 + CC(IDO,K,3)
521  106 CONTINUE
522  107 RETURN
523      END
524CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
525C
526C      File:           radf5.f
527C
528C      Library:        FFTPACK 4.1
529C
530C      Author:         Paul N. Swarztrauber
531C                      National Center for Atmospheric Research
532C                      PO 3000, Boulder, Colorado
533C
534C      Date:           Wed Mar 29 18:31:13 MST 1995
535C
536C      Description:    Lower-level auxiliary routine
537C
538      SUBROUTINE RADF5(IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
539      DOUBLE PRECISION TR11
540      DOUBLE PRECISION TI11
541      DOUBLE PRECISION TR12
542      DOUBLE PRECISION TI12
543      DOUBLE PRECISION CR2
544      DOUBLE PRECISION CI5
545      DOUBLE PRECISION CR3
546      DOUBLE PRECISION CI4
547      DOUBLE PRECISION DR2
548      DOUBLE PRECISION DI2
549      DOUBLE PRECISION DR3
550      DOUBLE PRECISION DI3
551      DOUBLE PRECISION DR4
552      DOUBLE PRECISION DI4
553      DOUBLE PRECISION DR5
554      DOUBLE PRECISION DI5
555      DOUBLE PRECISION CR5
556      DOUBLE PRECISION CI2
557      DOUBLE PRECISION CR4
558      DOUBLE PRECISION CI3
559      DOUBLE PRECISION TR2
560      DOUBLE PRECISION TI2
561      DOUBLE PRECISION TR3
562      DOUBLE PRECISION TI3
563      DOUBLE PRECISION TR5
564      DOUBLE PRECISION TI5
565      DOUBLE PRECISION TR4
566      DOUBLE PRECISION TI4
567      DOUBLE PRECISION CC(IDO,L1,5),CH(IDO,5,L1),WA1(*),WA2(*),WA3(*),
568     +                 WA4(*)
569      DATA TR11,TI11,TR12,TI12/.309016994374947D0,.951056516295154D0,
570     +     -.809016994374947D0,.587785252292473D0/
571C
572C FFTPACK 5.0 auxiliary routine
573C
574      DO 101 K = 1,L1
575          CR2 = CC(1,K,5) + CC(1,K,2)
576          CI5 = CC(1,K,5) - CC(1,K,2)
577          CR3 = CC(1,K,4) + CC(1,K,3)
578          CI4 = CC(1,K,4) - CC(1,K,3)
579          CH(1,1,K) = CC(1,K,1) + CR2 + CR3
580          CH(IDO,2,K) = CC(1,K,1) + TR11*CR2 + TR12*CR3
581          CH(1,3,K) = TI11*CI5 + TI12*CI4
582          CH(IDO,4,K) = CC(1,K,1) + TR12*CR2 + TR11*CR3
583          CH(1,5,K) = TI12*CI5 - TI11*CI4
584  101 CONTINUE
585      IF (IDO.EQ.1) RETURN
586      IDP2 = IDO + 2
587      DO 103 K = 1,L1
588          DO 102 I = 3,IDO,2
589              IC = IDP2 - I
590              DR2 = WA1(I-2)*CC(I-1,K,2) + WA1(I-1)*CC(I,K,2)
591              DI2 = WA1(I-2)*CC(I,K,2) - WA1(I-1)*CC(I-1,K,2)
592              DR3 = WA2(I-2)*CC(I-1,K,3) + WA2(I-1)*CC(I,K,3)
593              DI3 = WA2(I-2)*CC(I,K,3) - WA2(I-1)*CC(I-1,K,3)
594              DR4 = WA3(I-2)*CC(I-1,K,4) + WA3(I-1)*CC(I,K,4)
595              DI4 = WA3(I-2)*CC(I,K,4) - WA3(I-1)*CC(I-1,K,4)
596              DR5 = WA4(I-2)*CC(I-1,K,5) + WA4(I-1)*CC(I,K,5)
597              DI5 = WA4(I-2)*CC(I,K,5) - WA4(I-1)*CC(I-1,K,5)
598              CR2 = DR2 + DR5
599              CI5 = DR5 - DR2
600              CR5 = DI2 - DI5
601              CI2 = DI2 + DI5
602              CR3 = DR3 + DR4
603              CI4 = DR4 - DR3
604              CR4 = DI3 - DI4
605              CI3 = DI3 + DI4
606              CH(I-1,1,K) = CC(I-1,K,1) + CR2 + CR3
607              CH(I,1,K) = CC(I,K,1) + CI2 + CI3
608              TR2 = CC(I-1,K,1) + TR11*CR2 + TR12*CR3
609              TI2 = CC(I,K,1) + TR11*CI2 + TR12*CI3
610              TR3 = CC(I-1,K,1) + TR12*CR2 + TR11*CR3
611              TI3 = CC(I,K,1) + TR12*CI2 + TR11*CI3
612              TR5 = TI11*CR5 + TI12*CR4
613              TI5 = TI11*CI5 + TI12*CI4
614              TR4 = TI12*CR5 - TI11*CR4
615              TI4 = TI12*CI5 - TI11*CI4
616              CH(I-1,3,K) = TR2 + TR5
617              CH(IC-1,2,K) = TR2 - TR5
618              CH(I,3,K) = TI2 + TI5
619              CH(IC,2,K) = TI5 - TI2
620              CH(I-1,5,K) = TR3 + TR4
621              CH(IC-1,4,K) = TR3 - TR4
622              CH(I,5,K) = TI3 + TI4
623              CH(IC,4,K) = TI4 - TI3
624  102     CONTINUE
625  103 CONTINUE
626      RETURN
627      END
628CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
629C
630C      File:           radfg.f
631C
632C      Library:        FFTPACK 4.1
633C
634C      Author:         Paul N. Swarztrauber
635C                      National Center for Atmospheric Research
636C                      PO 3000, Boulder, Colorado
637C
638C      Date:           Wed Mar 29 18:31:13 MST 1995
639C
640C      Description:    Lower-level auxiliary routine
641C
642      SUBROUTINE RADFG(IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
643      DOUBLE PRECISION TPI
644      DOUBLE PRECISION PIMACH
645      DOUBLE PRECISION DUM
646      DOUBLE PRECISION ARG
647      DOUBLE PRECISION DCP
648      DOUBLE PRECISION DSP
649      DOUBLE PRECISION AR1
650      DOUBLE PRECISION AI1
651      DOUBLE PRECISION AR1H
652      DOUBLE PRECISION DC2
653      DOUBLE PRECISION DS2
654      DOUBLE PRECISION AR2
655      DOUBLE PRECISION AI2
656      DOUBLE PRECISION AR2H
657      DOUBLE PRECISION CH(IDO,L1,IP),CC(IDO,IP,L1),C1(IDO,L1,IP),
658     +                 C2(IDL1,IP),CH2(IDL1,IP),WA(*)
659C
660C FFTPACK 5.0 auxiliary routine
661C
662      TPI = 2.0D0*PIMACH(DUM)
663      ARG = TPI/DBLE(IP)
664      DCP = COS(ARG)
665      DSP = SIN(ARG)
666      IPPH = (IP+1)/2
667      IPP2 = IP + 2
668      IDP2 = IDO + 2
669      NBD = (IDO-1)/2
670      IF (IDO.EQ.1) GO TO 119
671      DO 101 IK = 1,IDL1
672          CH2(IK,1) = C2(IK,1)
673  101 CONTINUE
674      DO 103 J = 2,IP
675          DO 102 K = 1,L1
676              CH(1,K,J) = C1(1,K,J)
677  102     CONTINUE
678  103 CONTINUE
679      IF (NBD.GT.L1) GO TO 107
680      IS = -IDO
681      DO 106 J = 2,IP
682          IS = IS + IDO
683          IDIJ = IS
684          DO 105 I = 3,IDO,2
685              IDIJ = IDIJ + 2
686              DO 104 K = 1,L1
687                  CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J) +
688     +                          WA(IDIJ)*C1(I,K,J)
689                  CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J) -
690     +                        WA(IDIJ)*C1(I-1,K,J)
691  104         CONTINUE
692  105     CONTINUE
693  106 CONTINUE
694      GO TO 111
695  107 IS = -IDO
696      DO 110 J = 2,IP
697          IS = IS + IDO
698          DO 109 K = 1,L1
699              IDIJ = IS
700              DO 108 I = 3,IDO,2
701                  IDIJ = IDIJ + 2
702                  CH(I-1,K,J) = WA(IDIJ-1)*C1(I-1,K,J) +
703     +                          WA(IDIJ)*C1(I,K,J)
704                  CH(I,K,J) = WA(IDIJ-1)*C1(I,K,J) -
705     +                        WA(IDIJ)*C1(I-1,K,J)
706  108         CONTINUE
707  109     CONTINUE
708  110 CONTINUE
709  111 IF (NBD.LT.L1) GO TO 115
710      DO 114 J = 2,IPPH
711          JC = IPP2 - J
712          DO 113 K = 1,L1
713              DO 112 I = 3,IDO,2
714                  C1(I-1,K,J) = CH(I-1,K,J) + CH(I-1,K,JC)
715                  C1(I-1,K,JC) = CH(I,K,J) - CH(I,K,JC)
716                  C1(I,K,J) = CH(I,K,J) + CH(I,K,JC)
717                  C1(I,K,JC) = CH(I-1,K,JC) - CH(I-1,K,J)
718  112         CONTINUE
719  113     CONTINUE
720  114 CONTINUE
721      GO TO 121
722  115 DO 118 J = 2,IPPH
723          JC = IPP2 - J
724          DO 117 I = 3,IDO,2
725              DO 116 K = 1,L1
726                  C1(I-1,K,J) = CH(I-1,K,J) + CH(I-1,K,JC)
727                  C1(I-1,K,JC) = CH(I,K,J) - CH(I,K,JC)
728                  C1(I,K,J) = CH(I,K,J) + CH(I,K,JC)
729                  C1(I,K,JC) = CH(I-1,K,JC) - CH(I-1,K,J)
730  116         CONTINUE
731  117     CONTINUE
732  118 CONTINUE
733      GO TO 121
734  119 DO 120 IK = 1,IDL1
735          C2(IK,1) = CH2(IK,1)
736  120 CONTINUE
737  121 DO 123 J = 2,IPPH
738          JC = IPP2 - J
739          DO 122 K = 1,L1
740              C1(1,K,J) = CH(1,K,J) + CH(1,K,JC)
741              C1(1,K,JC) = CH(1,K,JC) - CH(1,K,J)
742  122     CONTINUE
743  123 CONTINUE
744C
745      AR1 = 1.D0
746      AI1 = 0.D0
747      DO 127 L = 2,IPPH
748          LC = IPP2 - L
749          AR1H = DCP*AR1 - DSP*AI1
750          AI1 = DCP*AI1 + DSP*AR1
751          AR1 = AR1H
752          DO 124 IK = 1,IDL1
753              CH2(IK,L) = C2(IK,1) + AR1*C2(IK,2)
754              CH2(IK,LC) = AI1*C2(IK,IP)
755  124     CONTINUE
756          DC2 = AR1
757          DS2 = AI1
758          AR2 = AR1
759          AI2 = AI1
760          DO 126 J = 3,IPPH
761              JC = IPP2 - J
762              AR2H = DC2*AR2 - DS2*AI2
763              AI2 = DC2*AI2 + DS2*AR2
764              AR2 = AR2H
765              DO 125 IK = 1,IDL1
766                  CH2(IK,L) = CH2(IK,L) + AR2*C2(IK,J)
767                  CH2(IK,LC) = CH2(IK,LC) + AI2*C2(IK,JC)
768  125         CONTINUE
769  126     CONTINUE
770  127 CONTINUE
771      DO 129 J = 2,IPPH
772          DO 128 IK = 1,IDL1
773              CH2(IK,1) = CH2(IK,1) + C2(IK,J)
774  128     CONTINUE
775  129 CONTINUE
776C
777      IF (IDO.LT.L1) GO TO 132
778      DO 131 K = 1,L1
779          DO 130 I = 1,IDO
780              CC(I,1,K) = CH(I,K,1)
781  130     CONTINUE
782  131 CONTINUE
783      GO TO 135
784  132 DO 134 I = 1,IDO
785          DO 133 K = 1,L1
786              CC(I,1,K) = CH(I,K,1)
787  133     CONTINUE
788  134 CONTINUE
789  135 DO 137 J = 2,IPPH
790          JC = IPP2 - J
791          J2 = J + J
792          DO 136 K = 1,L1
793              CC(IDO,J2-2,K) = CH(1,K,J)
794              CC(1,J2-1,K) = CH(1,K,JC)
795  136     CONTINUE
796  137 CONTINUE
797      IF (IDO.EQ.1) RETURN
798      IF (NBD.LT.L1) GO TO 141
799      DO 140 J = 2,IPPH
800          JC = IPP2 - J
801          J2 = J + J
802          DO 139 K = 1,L1
803              DO 138 I = 3,IDO,2
804                  IC = IDP2 - I
805                  CC(I-1,J2-1,K) = CH(I-1,K,J) + CH(I-1,K,JC)
806                  CC(IC-1,J2-2,K) = CH(I-1,K,J) - CH(I-1,K,JC)
807                  CC(I,J2-1,K) = CH(I,K,J) + CH(I,K,JC)
808                  CC(IC,J2-2,K) = CH(I,K,JC) - CH(I,K,J)
809  138         CONTINUE
810  139     CONTINUE
811  140 CONTINUE
812      RETURN
813  141 DO 144 J = 2,IPPH
814          JC = IPP2 - J
815          J2 = J + J
816          DO 143 I = 3,IDO,2
817              IC = IDP2 - I
818              DO 142 K = 1,L1
819                  CC(I-1,J2-1,K) = CH(I-1,K,J) + CH(I-1,K,JC)
820                  CC(IC-1,J2-2,K) = CH(I-1,K,J) - CH(I-1,K,JC)
821                  CC(I,J2-1,K) = CH(I,K,J) + CH(I,K,JC)
822                  CC(IC,J2-2,K) = CH(I,K,JC) - CH(I,K,J)
823  142         CONTINUE
824  143     CONTINUE
825  144 CONTINUE
826      RETURN
827      END
828