NCEPLIBS-sp  2.5.0
fftpack.F
Go to the documentation of this file.
1 C> @file
2 C> @brief A concatination of the (FFTPACK)[https://netlib.org/fftpack/] library code.
3 C>
4 C> FFTPACK is a package of Fortran subprograms for the fast Fourier
5 C> transform of periodic and other symmetric sequences. It includes
6 C> complex, real, sine, cosine, and quarter-wave transforms.
7 C>
8 C>Reference:
9 C>- P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
10 C>(G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.
11 C>
12 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
13 
14 C> dcrft
15 C>
16 C> @param init
17 C> @param x
18 C> @param ldx
19 C> @param y
20 C> @param ldy
21 C> @param n
22 C> @param m
23 C> @param isign
24 C> @param scale
25 C> @param table
26 C> @param n1
27 C> @param wrk
28 C> @param n2
29 C> @param z
30 C> @param nz
31 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
32  SUBROUTINE dcrft(init,x,ldx,y,ldy,n,m,isign,scale,
33  & table,n1,wrk,n2,z,nz)
34 
35  implicit none
36  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
37  real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
38 
39  IF (init.ne.0) THEN
40  CALL rffti(n,table)
41  ELSE
42 !OCL NOVREC
43  DO j=1,m
44  y(1,j)=x(1,j)
45  DO i=2,n
46  y(i,j)=x(i+1,j)
47  ENDDO
48  CALL rfftb(n,y(1,j),table)
49  DO i=1,n
50  y(i,j)=scale*y(i,j)
51  ENDDO
52  ENDDO
53  ENDIF
54 
55  RETURN
56  END
57 
58 C> scrft
59 C>
60 C> @param init
61 C> @param x
62 C> @param ldx
63 C> @param y
64 C> @param ldy
65 C> @param n
66 C> @param m
67 C> @param isign
68 C> @param scale
69 C> @param table
70 C> @param n1
71 C> @param wrk
72 C> @param n2
73 C> @param z
74 C> @param nz
75 C>
76 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
77 
78  SUBROUTINE scrft(init,x,ldx,y,ldy,n,m,isign,scale,
79  & table,n1,wrk,n2,z,nz)
80 
81  implicit none
82  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
83  real x(2*ldx,*),y(ldy,*),scale,table(44002),wrk,z
84 
85  IF (init.ne.0) THEN
86  CALL rffti(n,table)
87  ELSE
88 !OCL NOVREC
89  DO j=1,m
90  y(1,j)=x(1,j)
91  DO i=2,n
92  y(i,j)=x(i+1,j)
93  ENDDO
94  CALL rfftb(n,y(1,j),table)
95  DO i=1,n
96  y(i,j)=scale*y(i,j)
97  ENDDO
98  ENDDO
99  ENDIF
100 
101  RETURN
102  END
103 
104 C> csfft
105 C>
106 C> @param isign
107 C> @param n
108 C> @param scale
109 C> @param x
110 C> @param y
111 C> @param table
112 C> @param work
113 C> @param isys
114 C>
115 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
116  SUBROUTINE csfft(isign,n,scale,x,y,table,work,isys)
117 
118  implicit none
119  integer isign,n,isys,i
120  real scale,x(*),y(*),table(*),work(*)
121 
122  IF (isign.eq.0) THEN
123  CALL rffti(n,table)
124  ENDIF
125  IF (isign.eq.1) THEN
126  y(1)=x(1)
127  DO i=2,n
128  y(i)=x(i+1)
129  ENDDO
130  CALL rfftb(n,y,table)
131  DO i=1,n
132  y(i)=scale*y(i)
133  ENDDO
134  ENDIF
135 
136  RETURN
137  END
138 
139 C> drcft
140 C>
141 C> @param init
142 C> @param x
143 C> @param ldx
144 C> @param y
145 C> @param ldy
146 C> @param n
147 C> @param m
148 C> @param isign
149 C> @param scale
150 C> @param table
151 C> @param n1
152 C> @param wrk
153 C> @param n2
154 C> @param z
155 C> @param nz
156 C>
157 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
158  SUBROUTINE drcft(init,x,ldx,y,ldy,n,m,isign,scale,
159  & table,n1,wrk,n2,z,nz)
160 
161  implicit none
162  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
163  real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
164 
165  IF (init.ne.0) THEN
166  CALL rffti(n,table)
167  ELSE
168  DO j=1,m
169  DO i=1,n
170  y(i,j)=x(i,j)
171  ENDDO
172  CALL rfftf(n,y(1,j),table)
173  DO i=1,n
174  y(i,j)=scale*y(i,j)
175  ENDDO
176  DO i=n,2,-1
177  y(i+1,j)=y(i,j)
178  ENDDO
179  y(2,j)=0.
180 C 01/17/2013 vvvvvvvvvvvvv E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
181  y(n+2,j) = 0.
182  ENDDO
183  ENDIF
184 
185  RETURN
186  END
187 
188 C> srcft
189 C>
190 C> @param init
191 C> @param x
192 C> @param ldx
193 C> @param y
194 C> @param ldy
195 C> @param n
196 C> @param m
197 C> @param isign
198 C> @param scale
199 C> @param table
200 C> @param n1
201 C> @param wrk
202 C> @param n2
203 C> @param z
204 C> @param nz
205 C>
206 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
207  SUBROUTINE srcft(init,x,ldx,y,ldy,n,m,isign,scale,
208  & table,n1,wrk,n2,z,nz)
209 
210  implicit none
211  integer init,ldx,ldy,n,m,isign,n1,n2,nz,i,j
212  real x(ldx,*),y(2*ldy,*),scale,table(44002),wrk,z
213 
214  IF (init.ne.0) THEN
215  CALL rffti(n,table)
216  ELSE
217  DO j=1,m
218  DO i=1,n
219  y(i,j)=x(i,j)
220  ENDDO
221  CALL rfftf(n,y(1,j),table)
222  DO i=1,n
223  y(i,j)=scale*y(i,j)
224  ENDDO
225  DO i=n,2,-1
226  y(i+1,j)=y(i,j)
227  ENDDO
228  y(2,j)=0.
229  y(n+2,j) = 0.
230 C 01/17/2013 ^^^^^^^^^^E.Mirvis added ver 2.0.1 by S.Moorthi request. No +|- demo.
231  ENDDO
232  ENDIF
233 
234  RETURN
235  END
236 
237 C> scfft
238 C>
239 C> @param isign
240 C> @param n
241 C> @param scale
242 C> @param x
243 C> @param y
244 C> @param table
245 C> @param work
246 C> @param isys
247 C>
248 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
249  SUBROUTINE scfft(isign,n,scale,x,y,table,work,isys)
250 
251  implicit none
252  integer isign,n,isys,i
253  real scale,x(*),y(*),table(*),work(*)
254 
255  IF (isign.eq.0) THEN
256  CALL rffti(n,table)
257  ENDIF
258  IF (isign.eq.-1) THEN
259  DO i=1,n
260  y(i)=x(i)
261  ENDDO
262  CALL rfftf(n,y,table)
263  DO i=1,n
264  y(i)=scale*y(i)
265  ENDDO
266  DO i=n,2,-1
267  y(i+1)=y(i)
268  ENDDO
269  y(2)=0.
270  ENDIF
271 
272  RETURN
273  END
274 
275 C> RFFTF
276 C>
277 C> @param N
278 C> @param R
279 C> @param WSAVE
280 C>
281 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
282  SUBROUTINE rfftf (N,R,WSAVE)
283  dimension r(1) ,wsave(1)
284  IF (n .EQ. 1) RETURN
285  CALL rfftf1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
286  RETURN
287  END
288 
289 C> RFFTB
290 C>
291 C> @param N
292 C> @param R
293 C> @param WSAVE
294 C>
295 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
296  SUBROUTINE rfftb (N,R,WSAVE)
297  dimension r(1) ,wsave(1)
298  IF (n .EQ. 1) RETURN
299  CALL rfftb1 (n,r,wsave,wsave(n+1),wsave(2*n+1))
300  RETURN
301  END
302 
303 C> RFFTI
304 C>
305 C> @param N
306 C> @param WSAVE
307 C>
308 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
309  SUBROUTINE rffti (N,WSAVE)
310  dimension wsave(1)
311  IF (n .EQ. 1) RETURN
312  CALL rffti1 (n,wsave(n+1),wsave(2*n+1))
313  RETURN
314  END
315 
316 C> RFFTB1
317 C>
318 C> @param N
319 C> @param C
320 C> @param CH
321 C> @param WA
322 C> @param IFAC
323 C>
324 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
325  SUBROUTINE rfftb1 (N,C,CH,WA,IFAC)
326  dimension ch(1) ,c(1) ,wa(1) ,ifac(*)
327  nf = ifac(2)
328  na = 0
329  l1 = 1
330  iw = 1
331  DO 116 k1=1,nf
332  ip = ifac(k1+2)
333  l2 = ip*l1
334  ido = n/l2
335  idl1 = ido*l1
336  IF (ip .NE. 4) GO TO 103
337  ix2 = iw+ido
338  ix3 = ix2+ido
339  IF (na .NE. 0) GO TO 101
340  CALL radb4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
341  GO TO 102
342  101 CALL radb4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
343  102 na = 1-na
344  GO TO 115
345  103 IF (ip .NE. 2) GO TO 106
346  IF (na .NE. 0) GO TO 104
347  CALL radb2 (ido,l1,c,ch,wa(iw))
348  GO TO 105
349  104 CALL radb2 (ido,l1,ch,c,wa(iw))
350  105 na = 1-na
351  GO TO 115
352  106 IF (ip .NE. 3) GO TO 109
353  ix2 = iw+ido
354  IF (na .NE. 0) GO TO 107
355  CALL radb3 (ido,l1,c,ch,wa(iw),wa(ix2))
356  GO TO 108
357  107 CALL radb3 (ido,l1,ch,c,wa(iw),wa(ix2))
358  108 na = 1-na
359  GO TO 115
360  109 IF (ip .NE. 5) GO TO 112
361  ix2 = iw+ido
362  ix3 = ix2+ido
363  ix4 = ix3+ido
364  IF (na .NE. 0) GO TO 110
365  CALL radb5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
366  GO TO 111
367  110 CALL radb5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
368  111 na = 1-na
369  GO TO 115
370  112 IF (na .NE. 0) GO TO 113
371  CALL radbg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
372  GO TO 114
373  113 CALL radbg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
374  114 IF (ido .EQ. 1) na = 1-na
375  115 l1 = l2
376  iw = iw+(ip-1)*ido
377  116 CONTINUE
378  IF (na .EQ. 0) RETURN
379  DO 117 i=1,n
380  c(i) = ch(i)
381  117 CONTINUE
382  RETURN
383  END
384 
385 C> RFFTF1
386 C>
387 C> @param N
388 C> @param C
389 C> @param CH
390 C> @param WA
391 C> @param IFAC
392 C>
393 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
394  SUBROUTINE rfftf1 (N,C,CH,WA,IFAC)
395  dimension ch(1) ,c(1) ,wa(1) ,ifac(*)
396  nf = ifac(2)
397  na = 1
398  l2 = n
399  iw = n
400  DO 111 k1=1,nf
401  kh = nf-k1
402  ip = ifac(kh+3)
403  l1 = l2/ip
404  ido = n/l2
405  idl1 = ido*l1
406  iw = iw-(ip-1)*ido
407  na = 1-na
408  IF (ip .NE. 4) GO TO 102
409  ix2 = iw+ido
410  ix3 = ix2+ido
411  IF (na .NE. 0) GO TO 101
412  CALL radf4 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3))
413  GO TO 110
414  101 CALL radf4 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3))
415  GO TO 110
416  102 IF (ip .NE. 2) GO TO 104
417  IF (na .NE. 0) GO TO 103
418  CALL radf2 (ido,l1,c,ch,wa(iw))
419  GO TO 110
420  103 CALL radf2 (ido,l1,ch,c,wa(iw))
421  GO TO 110
422  104 IF (ip .NE. 3) GO TO 106
423  ix2 = iw+ido
424  IF (na .NE. 0) GO TO 105
425  CALL radf3 (ido,l1,c,ch,wa(iw),wa(ix2))
426  GO TO 110
427  105 CALL radf3 (ido,l1,ch,c,wa(iw),wa(ix2))
428  GO TO 110
429  106 IF (ip .NE. 5) GO TO 108
430  ix2 = iw+ido
431  ix3 = ix2+ido
432  ix4 = ix3+ido
433  IF (na .NE. 0) GO TO 107
434  CALL radf5 (ido,l1,c,ch,wa(iw),wa(ix2),wa(ix3),wa(ix4))
435  GO TO 110
436  107 CALL radf5 (ido,l1,ch,c,wa(iw),wa(ix2),wa(ix3),wa(ix4))
437  GO TO 110
438  108 IF (ido .EQ. 1) na = 1-na
439  IF (na .NE. 0) GO TO 109
440  CALL radfg (ido,ip,l1,idl1,c,c,c,ch,ch,wa(iw))
441  na = 1
442  GO TO 110
443  109 CALL radfg (ido,ip,l1,idl1,ch,ch,ch,c,c,wa(iw))
444  na = 0
445  110 l2 = l1
446  111 CONTINUE
447  IF (na .EQ. 1) RETURN
448  DO 112 i=1,n
449  c(i) = ch(i)
450  112 CONTINUE
451  RETURN
452  END
453 
454 C> RFFTI1
455 C>
456 C> @param N
457 C> @param WA
458 C> @param IFAC
459 C>
460 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
461  SUBROUTINE rffti1 (N,WA,IFAC)
462  dimension wa(1) ,ifac(*) ,ntryh(4)
463  DATA ntryh(1),ntryh(2),ntryh(3),ntryh(4)/4,2,3,5/
464  nl = n
465  nf = 0
466  j = 0
467  101 j = j+1
468  IF (j-4) 102,102,103
469  102 ntry = ntryh(j)
470  GO TO 104
471  103 ntry = ntry+2
472  104 nq = nl/ntry
473  nr = nl-ntry*nq
474  IF (nr) 101,105,101
475  105 nf = nf+1
476  ifac(nf+2) = ntry
477  nl = nq
478  IF (ntry .NE. 2) GO TO 107
479  IF (nf .EQ. 1) GO TO 107
480  DO 106 i=2,nf
481  ib = nf-i+2
482  ifac(ib+2) = ifac(ib+1)
483  106 CONTINUE
484  ifac(3) = 2
485  107 IF (nl .NE. 1) GO TO 104
486  ifac(1) = n
487  ifac(2) = nf
488  tpi = 6.28318530717959
489  argh = tpi/float(n)
490  is = 0
491  nfm1 = nf-1
492  l1 = 1
493  IF (nfm1 .EQ. 0) RETURN
494 !OCL NOVREC
495  DO 110 k1=1,nfm1
496  ip = ifac(k1+2)
497  ld = 0
498  l2 = l1*ip
499  ido = n/l2
500  ipm = ip-1
501  DO 109 j=1,ipm
502  ld = ld+l1
503  i = is
504  argld = float(ld)*argh
505  fi = 0
506 !OCL SCALAR
507  DO 108 ii=3,ido,2
508  i = i+2
509  fi = fi+1
510  arg = fi*argld
511  wa(i-1) = cos(arg)
512  wa(i) = sin(arg)
513  108 CONTINUE
514  is = is+ido
515  109 CONTINUE
516  l1 = l2
517  110 CONTINUE
518  RETURN
519  END
520 
521 C> RADB2
522 C>
523 C> @param IDO
524 C> @param L1
525 C> @param CC
526 C> @param CH
527 C> @param WA1
528 C>
529 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
530  SUBROUTINE radb2 (IDO,L1,CC,CH,WA1)
531  dimension cc(ido,2,l1) ,ch(ido,l1,2) ,
532  1 wa1(1)
533  DO 101 k=1,l1
534  ch(1,k,1) = cc(1,1,k)+cc(ido,2,k)
535  ch(1,k,2) = cc(1,1,k)-cc(ido,2,k)
536  101 CONTINUE
537  IF (ido-2) 107,105,102
538  102 idp2 = ido+2
539 !OCL NOVREC
540  DO 104 k=1,l1
541  DO 103 i=3,ido,2
542  ic = idp2-i
543  ch(i-1,k,1) = cc(i-1,1,k)+cc(ic-1,2,k)
544  tr2 = cc(i-1,1,k)-cc(ic-1,2,k)
545  ch(i,k,1) = cc(i,1,k)-cc(ic,2,k)
546  ti2 = cc(i,1,k)+cc(ic,2,k)
547  ch(i-1,k,2) = wa1(i-2)*tr2-wa1(i-1)*ti2
548  ch(i,k,2) = wa1(i-2)*ti2+wa1(i-1)*tr2
549  103 CONTINUE
550  104 CONTINUE
551  IF (mod(ido,2) .EQ. 1) RETURN
552  105 DO 106 k=1,l1
553  ch(ido,k,1) = cc(ido,1,k)+cc(ido,1,k)
554  ch(ido,k,2) = -(cc(1,2,k)+cc(1,2,k))
555  106 CONTINUE
556  107 RETURN
557  END
558 
559 C> RADB3
560 C>
561 C> @param IDO
562 C> @param L1
563 C> @param CC
564 C> @param CH
565 C> @param WA1
566 C> @param WA2
567 C>
568 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
569  SUBROUTINE radb3 (IDO,L1,CC,CH,WA1,WA2)
570  dimension cc(ido,3,l1) ,ch(ido,l1,3) ,
571  1 wa1(1) ,wa2(1)
572  DATA taur,taui /-.5,.866025403784439/
573  DO 101 k=1,l1
574  tr2 = cc(ido,2,k)+cc(ido,2,k)
575  cr2 = cc(1,1,k)+taur*tr2
576  ch(1,k,1) = cc(1,1,k)+tr2
577  ci3 = taui*(cc(1,3,k)+cc(1,3,k))
578  ch(1,k,2) = cr2-ci3
579  ch(1,k,3) = cr2+ci3
580  101 CONTINUE
581  IF (ido .EQ. 1) RETURN
582  idp2 = ido+2
583 !OCL NOVREC
584  DO 103 k=1,l1
585  DO 102 i=3,ido,2
586  ic = idp2-i
587  tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
588  cr2 = cc(i-1,1,k)+taur*tr2
589  ch(i-1,k,1) = cc(i-1,1,k)+tr2
590  ti2 = cc(i,3,k)-cc(ic,2,k)
591  ci2 = cc(i,1,k)+taur*ti2
592  ch(i,k,1) = cc(i,1,k)+ti2
593  cr3 = taui*(cc(i-1,3,k)-cc(ic-1,2,k))
594  ci3 = taui*(cc(i,3,k)+cc(ic,2,k))
595  dr2 = cr2-ci3
596  dr3 = cr2+ci3
597  di2 = ci2+cr3
598  di3 = ci2-cr3
599  ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
600  ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
601  ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
602  ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
603  102 CONTINUE
604  103 CONTINUE
605  RETURN
606  END
607 
608 C> RADB4
609 C>
610 C> @param IDO
611 C> @param L1
612 C> @param CC
613 C> @param CH
614 C> @param WA1
615 C> @param WA2
616 C> @param WA3
617 C>
618 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
619  SUBROUTINE radb4 (IDO,L1,CC,CH,WA1,WA2,WA3)
620  dimension cc(ido,4,l1) ,ch(ido,l1,4) ,
621  1 wa1(1) ,wa2(1) ,wa3(1)
622  DATA sqrt2 /1.414213562373095/
623  DO 101 k=1,l1
624  tr1 = cc(1,1,k)-cc(ido,4,k)
625  tr2 = cc(1,1,k)+cc(ido,4,k)
626  tr3 = cc(ido,2,k)+cc(ido,2,k)
627  tr4 = cc(1,3,k)+cc(1,3,k)
628  ch(1,k,1) = tr2+tr3
629  ch(1,k,2) = tr1-tr4
630  ch(1,k,3) = tr2-tr3
631  ch(1,k,4) = tr1+tr4
632  101 CONTINUE
633  IF (ido-2) 107,105,102
634  102 idp2 = ido+2
635 !OCL NOVREC
636  DO 104 k=1,l1
637  DO 103 i=3,ido,2
638  ic = idp2-i
639  ti1 = cc(i,1,k)+cc(ic,4,k)
640  ti2 = cc(i,1,k)-cc(ic,4,k)
641  ti3 = cc(i,3,k)-cc(ic,2,k)
642  tr4 = cc(i,3,k)+cc(ic,2,k)
643  tr1 = cc(i-1,1,k)-cc(ic-1,4,k)
644  tr2 = cc(i-1,1,k)+cc(ic-1,4,k)
645  ti4 = cc(i-1,3,k)-cc(ic-1,2,k)
646  tr3 = cc(i-1,3,k)+cc(ic-1,2,k)
647  ch(i-1,k,1) = tr2+tr3
648  cr3 = tr2-tr3
649  ch(i,k,1) = ti2+ti3
650  ci3 = ti2-ti3
651  cr2 = tr1-tr4
652  cr4 = tr1+tr4
653  ci2 = ti1+ti4
654  ci4 = ti1-ti4
655  ch(i-1,k,2) = wa1(i-2)*cr2-wa1(i-1)*ci2
656  ch(i,k,2) = wa1(i-2)*ci2+wa1(i-1)*cr2
657  ch(i-1,k,3) = wa2(i-2)*cr3-wa2(i-1)*ci3
658  ch(i,k,3) = wa2(i-2)*ci3+wa2(i-1)*cr3
659  ch(i-1,k,4) = wa3(i-2)*cr4-wa3(i-1)*ci4
660  ch(i,k,4) = wa3(i-2)*ci4+wa3(i-1)*cr4
661  103 CONTINUE
662  104 CONTINUE
663  IF (mod(ido,2) .EQ. 1) RETURN
664  105 CONTINUE
665  DO 106 k=1,l1
666  ti1 = cc(1,2,k)+cc(1,4,k)
667  ti2 = cc(1,4,k)-cc(1,2,k)
668  tr1 = cc(ido,1,k)-cc(ido,3,k)
669  tr2 = cc(ido,1,k)+cc(ido,3,k)
670  ch(ido,k,1) = tr2+tr2
671  ch(ido,k,2) = sqrt2*(tr1-ti1)
672  ch(ido,k,3) = ti2+ti2
673  ch(ido,k,4) = -sqrt2*(tr1+ti1)
674  106 CONTINUE
675  107 RETURN
676  END
677 
678 C> RADB5
679 C>
680 C> @param IDO
681 C> @param L1
682 C> @param CC
683 C> @param CH
684 C> @param WA1
685 C> @param WA2
686 C> @param WA3
687 C> @param WA4
688 C>
689 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
690  SUBROUTINE radb5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
691  dimension cc(ido,5,l1) ,ch(ido,l1,5) ,
692  1 wa1(1) ,wa2(1) ,wa3(1) ,wa4(1)
693  DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
694  1-.809016994374947,.587785252292473/
695  DO 101 k=1,l1
696  ti5 = cc(1,3,k)+cc(1,3,k)
697  ti4 = cc(1,5,k)+cc(1,5,k)
698  tr2 = cc(ido,2,k)+cc(ido,2,k)
699  tr3 = cc(ido,4,k)+cc(ido,4,k)
700  ch(1,k,1) = cc(1,1,k)+tr2+tr3
701  cr2 = cc(1,1,k)+tr11*tr2+tr12*tr3
702  cr3 = cc(1,1,k)+tr12*tr2+tr11*tr3
703  ci5 = ti11*ti5+ti12*ti4
704  ci4 = ti12*ti5-ti11*ti4
705  ch(1,k,2) = cr2-ci5
706  ch(1,k,3) = cr3-ci4
707  ch(1,k,4) = cr3+ci4
708  ch(1,k,5) = cr2+ci5
709  101 CONTINUE
710  IF (ido .EQ. 1) RETURN
711  idp2 = ido+2
712  DO 103 k=1,l1
713  DO 102 i=3,ido,2
714  ic = idp2-i
715  ti5 = cc(i,3,k)+cc(ic,2,k)
716  ti2 = cc(i,3,k)-cc(ic,2,k)
717  ti4 = cc(i,5,k)+cc(ic,4,k)
718  ti3 = cc(i,5,k)-cc(ic,4,k)
719  tr5 = cc(i-1,3,k)-cc(ic-1,2,k)
720  tr2 = cc(i-1,3,k)+cc(ic-1,2,k)
721  tr4 = cc(i-1,5,k)-cc(ic-1,4,k)
722  tr3 = cc(i-1,5,k)+cc(ic-1,4,k)
723  ch(i-1,k,1) = cc(i-1,1,k)+tr2+tr3
724  ch(i,k,1) = cc(i,1,k)+ti2+ti3
725  cr2 = cc(i-1,1,k)+tr11*tr2+tr12*tr3
726  ci2 = cc(i,1,k)+tr11*ti2+tr12*ti3
727  cr3 = cc(i-1,1,k)+tr12*tr2+tr11*tr3
728  ci3 = cc(i,1,k)+tr12*ti2+tr11*ti3
729  cr5 = ti11*tr5+ti12*tr4
730  ci5 = ti11*ti5+ti12*ti4
731  cr4 = ti12*tr5-ti11*tr4
732  ci4 = ti12*ti5-ti11*ti4
733  dr3 = cr3-ci4
734  dr4 = cr3+ci4
735  di3 = ci3+cr4
736  di4 = ci3-cr4
737  dr5 = cr2+ci5
738  dr2 = cr2-ci5
739  di5 = ci2-cr5
740  di2 = ci2+cr5
741  ch(i-1,k,2) = wa1(i-2)*dr2-wa1(i-1)*di2
742  ch(i,k,2) = wa1(i-2)*di2+wa1(i-1)*dr2
743  ch(i-1,k,3) = wa2(i-2)*dr3-wa2(i-1)*di3
744  ch(i,k,3) = wa2(i-2)*di3+wa2(i-1)*dr3
745  ch(i-1,k,4) = wa3(i-2)*dr4-wa3(i-1)*di4
746  ch(i,k,4) = wa3(i-2)*di4+wa3(i-1)*dr4
747  ch(i-1,k,5) = wa4(i-2)*dr5-wa4(i-1)*di5
748  ch(i,k,5) = wa4(i-2)*di5+wa4(i-1)*dr5
749  102 CONTINUE
750  103 CONTINUE
751  RETURN
752  END
753 
754 C> RADBG
755 C>
756 C> @param IDO
757 C> @param IP
758 C> @param L1
759 C> @param IDL1
760 C> @param CC
761 C> @param C1
762 C> @param C2
763 C> @param CH
764 C> @param CH2
765 C> @param WA
766 C>
767 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
768  SUBROUTINE radbg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
769  dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
770  1 c1(ido,l1,ip) ,c2(idl1,ip),
771  2 ch2(idl1,ip) ,wa(1)
772  DATA tpi/6.28318530717959/
773  arg = tpi/float(ip)
774  dcp = cos(arg)
775  dsp = sin(arg)
776  idp2 = ido+2
777  nbd = (ido-1)/2
778  ipp2 = ip+2
779  ipph = (ip+1)/2
780  IF (ido .LT. l1) GO TO 103
781  DO 102 k=1,l1
782  DO 101 i=1,ido
783  ch(i,k,1) = cc(i,1,k)
784  101 CONTINUE
785  102 CONTINUE
786  GO TO 106
787  103 DO 105 i=1,ido
788  DO 104 k=1,l1
789  ch(i,k,1) = cc(i,1,k)
790  104 CONTINUE
791  105 CONTINUE
792 !OCL NOVREC
793  106 DO 108 j=2,ipph
794  jc = ipp2-j
795  j2 = j+j
796  DO 107 k=1,l1
797  ch(1,k,j) = cc(ido,j2-2,k)+cc(ido,j2-2,k)
798  ch(1,k,jc) = cc(1,j2-1,k)+cc(1,j2-1,k)
799  107 CONTINUE
800  108 CONTINUE
801  IF (ido .EQ. 1) GO TO 116
802  IF (nbd .LT. l1) GO TO 112
803 !OCL NOVREC
804  DO 111 j=2,ipph
805  jc = ipp2-j
806  DO 110 k=1,l1
807  DO 109 i=3,ido,2
808  ic = idp2-i
809  ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
810  ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
811  ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
812  ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
813  109 CONTINUE
814  110 CONTINUE
815  111 CONTINUE
816  GO TO 116
817  112 DO 115 j=2,ipph
818  jc = ipp2-j
819  DO 114 i=3,ido,2
820  ic = idp2-i
821  DO 113 k=1,l1
822  ch(i-1,k,j) = cc(i-1,2*j-1,k)+cc(ic-1,2*j-2,k)
823  ch(i-1,k,jc) = cc(i-1,2*j-1,k)-cc(ic-1,2*j-2,k)
824  ch(i,k,j) = cc(i,2*j-1,k)-cc(ic,2*j-2,k)
825  ch(i,k,jc) = cc(i,2*j-1,k)+cc(ic,2*j-2,k)
826  113 CONTINUE
827  114 CONTINUE
828  115 CONTINUE
829  116 ar1 = 1.
830  ai1 = 0.
831 !OCL NOVREC
832  DO 120 l=2,ipph
833  lc = ipp2-l
834  ar1h = dcp*ar1-dsp*ai1
835  ai1 = dcp*ai1+dsp*ar1
836  ar1 = ar1h
837  DO 117 ik=1,idl1
838  c2(ik,l) = ch2(ik,1)+ar1*ch2(ik,2)
839  c2(ik,lc) = ai1*ch2(ik,ip)
840  117 CONTINUE
841  dc2 = ar1
842  ds2 = ai1
843  ar2 = ar1
844  ai2 = ai1
845 !OCL NOVREC
846  DO 119 j=3,ipph
847  jc = ipp2-j
848  ar2h = dc2*ar2-ds2*ai2
849  ai2 = dc2*ai2+ds2*ar2
850  ar2 = ar2h
851  DO 118 ik=1,idl1
852  c2(ik,l) = c2(ik,l)+ar2*ch2(ik,j)
853  c2(ik,lc) = c2(ik,lc)+ai2*ch2(ik,jc)
854  118 CONTINUE
855  119 CONTINUE
856  120 CONTINUE
857 !OCL NOVREC
858  DO 122 j=2,ipph
859  DO 121 ik=1,idl1
860  ch2(ik,1) = ch2(ik,1)+ch2(ik,j)
861  121 CONTINUE
862  122 CONTINUE
863 !OCL NOVREC
864  DO 124 j=2,ipph
865  jc = ipp2-j
866  DO 123 k=1,l1
867  ch(1,k,j) = c1(1,k,j)-c1(1,k,jc)
868  ch(1,k,jc) = c1(1,k,j)+c1(1,k,jc)
869  123 CONTINUE
870  124 CONTINUE
871  IF (ido .EQ. 1) GO TO 132
872  IF (nbd .LT. l1) GO TO 128
873 !OCL NOVREC
874  DO 127 j=2,ipph
875  jc = ipp2-j
876  DO 126 k=1,l1
877  DO 125 i=3,ido,2
878  ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
879  ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
880  ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
881  ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
882  125 CONTINUE
883  126 CONTINUE
884  127 CONTINUE
885  GO TO 132
886  128 DO 131 j=2,ipph
887  jc = ipp2-j
888  DO 130 i=3,ido,2
889  DO 129 k=1,l1
890  ch(i-1,k,j) = c1(i-1,k,j)-c1(i,k,jc)
891  ch(i-1,k,jc) = c1(i-1,k,j)+c1(i,k,jc)
892  ch(i,k,j) = c1(i,k,j)+c1(i-1,k,jc)
893  ch(i,k,jc) = c1(i,k,j)-c1(i-1,k,jc)
894  129 CONTINUE
895  130 CONTINUE
896  131 CONTINUE
897  132 CONTINUE
898  IF (ido .EQ. 1) RETURN
899  DO 133 ik=1,idl1
900  c2(ik,1) = ch2(ik,1)
901  133 CONTINUE
902  DO 135 j=2,ip
903  DO 134 k=1,l1
904  c1(1,k,j) = ch(1,k,j)
905  134 CONTINUE
906  135 CONTINUE
907  IF (nbd .GT. l1) GO TO 139
908  is = -ido
909  DO 138 j=2,ip
910  is = is+ido
911  idij = is
912  DO 137 i=3,ido,2
913  idij = idij+2
914  DO 136 k=1,l1
915  c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
916  c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
917  136 CONTINUE
918  137 CONTINUE
919  138 CONTINUE
920  GO TO 143
921  139 is = -ido
922 !OCL NOVREC
923  DO 142 j=2,ip
924  is = is+ido
925  DO 141 k=1,l1
926  idij = is
927  DO 140 i=3,ido,2
928  idij = idij+2
929  c1(i-1,k,j) = wa(idij-1)*ch(i-1,k,j)-wa(idij)*ch(i,k,j)
930  c1(i,k,j) = wa(idij-1)*ch(i,k,j)+wa(idij)*ch(i-1,k,j)
931  140 CONTINUE
932  141 CONTINUE
933  142 CONTINUE
934  143 RETURN
935  END
936 
937 C> RADBG
938 C>
939 C> @param IDO
940 C> @param L1
941 C> @param CC
942 C> @param CH
943 C> @param WA1
944 C>
945 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
946  SUBROUTINE radf2 (IDO,L1,CC,CH,WA1)
947  dimension ch(ido,2,l1) ,cc(ido,l1,2) ,
948  1 wa1(1)
949  DO 101 k=1,l1
950  ch(1,1,k) = cc(1,k,1)+cc(1,k,2)
951  ch(ido,2,k) = cc(1,k,1)-cc(1,k,2)
952  101 CONTINUE
953  IF (ido-2) 107,105,102
954  102 idp2 = ido+2
955  DO 104 k=1,l1
956  DO 103 i=3,ido,2
957  ic = idp2-i
958  tr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
959  ti2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
960  ch(i,1,k) = cc(i,k,1)+ti2
961  ch(ic,2,k) = ti2-cc(i,k,1)
962  ch(i-1,1,k) = cc(i-1,k,1)+tr2
963  ch(ic-1,2,k) = cc(i-1,k,1)-tr2
964  103 CONTINUE
965  104 CONTINUE
966  IF (mod(ido,2) .EQ. 1) RETURN
967  105 DO 106 k=1,l1
968  ch(1,2,k) = -cc(ido,k,2)
969  ch(ido,1,k) = cc(ido,k,1)
970  106 CONTINUE
971  107 RETURN
972  END
973 
974 C> RADF3
975 C>
976 C> @param IDO
977 C> @param L1
978 C> @param CC
979 C> @param CH
980 C> @param WA1
981 C> @param WA2
982 C>
983 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
984  SUBROUTINE radf3 (IDO,L1,CC,CH,WA1,WA2)
985  dimension ch(ido,3,l1) ,cc(ido,l1,3) ,
986  1 wa1(1) ,wa2(1)
987  DATA taur,taui /-.5,.866025403784439/
988  DO 101 k=1,l1
989  cr2 = cc(1,k,2)+cc(1,k,3)
990  ch(1,1,k) = cc(1,k,1)+cr2
991  ch(1,3,k) = taui*(cc(1,k,3)-cc(1,k,2))
992  ch(ido,2,k) = cc(1,k,1)+taur*cr2
993  101 CONTINUE
994  IF (ido .EQ. 1) RETURN
995  idp2 = ido+2
996  DO 103 k=1,l1
997  DO 102 i=3,ido,2
998  ic = idp2-i
999  dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1000  di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1001  dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1002  di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1003  cr2 = dr2+dr3
1004  ci2 = di2+di3
1005  ch(i-1,1,k) = cc(i-1,k,1)+cr2
1006  ch(i,1,k) = cc(i,k,1)+ci2
1007  tr2 = cc(i-1,k,1)+taur*cr2
1008  ti2 = cc(i,k,1)+taur*ci2
1009  tr3 = taui*(di2-di3)
1010  ti3 = taui*(dr3-dr2)
1011  ch(i-1,3,k) = tr2+tr3
1012  ch(ic-1,2,k) = tr2-tr3
1013  ch(i,3,k) = ti2+ti3
1014  ch(ic,2,k) = ti3-ti2
1015  102 CONTINUE
1016  103 CONTINUE
1017  RETURN
1018  END
1019 
1020 C> RADF4
1021 C>
1022 C> @param IDO
1023 C> @param L1
1024 C> @param CC
1025 C> @param CH
1026 C> @param WA1
1027 C> @param WA2
1028 C> @param WA3
1029 C>
1030 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1031  SUBROUTINE radf4 (IDO,L1,CC,CH,WA1,WA2,WA3)
1032  dimension cc(ido,l1,4) ,ch(ido,4,l1) ,
1033  1 wa1(1) ,wa2(1) ,wa3(1)
1034  DATA hsqt2 /.7071067811865475/
1035  DO 101 k=1,l1
1036  tr1 = cc(1,k,2)+cc(1,k,4)
1037  tr2 = cc(1,k,1)+cc(1,k,3)
1038  ch(1,1,k) = tr1+tr2
1039  ch(ido,4,k) = tr2-tr1
1040  ch(ido,2,k) = cc(1,k,1)-cc(1,k,3)
1041  ch(1,3,k) = cc(1,k,4)-cc(1,k,2)
1042  101 CONTINUE
1043  IF (ido-2) 107,105,102
1044  102 idp2 = ido+2
1045 !OCL NOVREC
1046  DO 104 k=1,l1
1047  DO 103 i=3,ido,2
1048  ic = idp2-i
1049  cr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1050  ci2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1051  cr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1052  ci3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1053  cr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
1054  ci4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
1055  tr1 = cr2+cr4
1056  tr4 = cr4-cr2
1057  ti1 = ci2+ci4
1058  ti4 = ci2-ci4
1059  ti2 = cc(i,k,1)+ci3
1060  ti3 = cc(i,k,1)-ci3
1061  tr2 = cc(i-1,k,1)+cr3
1062  tr3 = cc(i-1,k,1)-cr3
1063  ch(i-1,1,k) = tr1+tr2
1064  ch(ic-1,4,k) = tr2-tr1
1065  ch(i,1,k) = ti1+ti2
1066  ch(ic,4,k) = ti1-ti2
1067  ch(i-1,3,k) = ti4+tr3
1068  ch(ic-1,2,k) = tr3-ti4
1069  ch(i,3,k) = tr4+ti3
1070  ch(ic,2,k) = tr4-ti3
1071  103 CONTINUE
1072  104 CONTINUE
1073  IF (mod(ido,2) .EQ. 1) RETURN
1074  105 CONTINUE
1075  DO 106 k=1,l1
1076  ti1 = -hsqt2*(cc(ido,k,2)+cc(ido,k,4))
1077  tr1 = hsqt2*(cc(ido,k,2)-cc(ido,k,4))
1078  ch(ido,1,k) = tr1+cc(ido,k,1)
1079  ch(ido,3,k) = cc(ido,k,1)-tr1
1080  ch(1,2,k) = ti1-cc(ido,k,3)
1081  ch(1,4,k) = ti1+cc(ido,k,3)
1082  106 CONTINUE
1083  107 RETURN
1084  END
1085 
1086 C> RADF5
1087 C>
1088 C> @param IDO
1089 C> @param L1
1090 C> @param CC
1091 C> @param CH
1092 C> @param WA1
1093 C> @param WA2
1094 C> @param WA3
1095 C> @param WA4
1096 C>
1097 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1098  SUBROUTINE radf5 (IDO,L1,CC,CH,WA1,WA2,WA3,WA4)
1099  dimension cc(ido,l1,5) ,ch(ido,5,l1) ,
1100  1 wa1(1) ,wa2(1) ,wa3(1) ,wa4(1)
1101  DATA tr11,ti11,tr12,ti12 /.309016994374947,.951056516295154,
1102  1-.809016994374947,.587785252292473/
1103  DO 101 k=1,l1
1104  cr2 = cc(1,k,5)+cc(1,k,2)
1105  ci5 = cc(1,k,5)-cc(1,k,2)
1106  cr3 = cc(1,k,4)+cc(1,k,3)
1107  ci4 = cc(1,k,4)-cc(1,k,3)
1108  ch(1,1,k) = cc(1,k,1)+cr2+cr3
1109  ch(ido,2,k) = cc(1,k,1)+tr11*cr2+tr12*cr3
1110  ch(1,3,k) = ti11*ci5+ti12*ci4
1111  ch(ido,4,k) = cc(1,k,1)+tr12*cr2+tr11*cr3
1112  ch(1,5,k) = ti12*ci5-ti11*ci4
1113  101 CONTINUE
1114  IF (ido .EQ. 1) RETURN
1115  idp2 = ido+2
1116  DO 103 k=1,l1
1117  DO 102 i=3,ido,2
1118  ic = idp2-i
1119  dr2 = wa1(i-2)*cc(i-1,k,2)+wa1(i-1)*cc(i,k,2)
1120  di2 = wa1(i-2)*cc(i,k,2)-wa1(i-1)*cc(i-1,k,2)
1121  dr3 = wa2(i-2)*cc(i-1,k,3)+wa2(i-1)*cc(i,k,3)
1122  di3 = wa2(i-2)*cc(i,k,3)-wa2(i-1)*cc(i-1,k,3)
1123  dr4 = wa3(i-2)*cc(i-1,k,4)+wa3(i-1)*cc(i,k,4)
1124  di4 = wa3(i-2)*cc(i,k,4)-wa3(i-1)*cc(i-1,k,4)
1125  dr5 = wa4(i-2)*cc(i-1,k,5)+wa4(i-1)*cc(i,k,5)
1126  di5 = wa4(i-2)*cc(i,k,5)-wa4(i-1)*cc(i-1,k,5)
1127  cr2 = dr2+dr5
1128  ci5 = dr5-dr2
1129  cr5 = di2-di5
1130  ci2 = di2+di5
1131  cr3 = dr3+dr4
1132  ci4 = dr4-dr3
1133  cr4 = di3-di4
1134  ci3 = di3+di4
1135  ch(i-1,1,k) = cc(i-1,k,1)+cr2+cr3
1136  ch(i,1,k) = cc(i,k,1)+ci2+ci3
1137  tr2 = cc(i-1,k,1)+tr11*cr2+tr12*cr3
1138  ti2 = cc(i,k,1)+tr11*ci2+tr12*ci3
1139  tr3 = cc(i-1,k,1)+tr12*cr2+tr11*cr3
1140  ti3 = cc(i,k,1)+tr12*ci2+tr11*ci3
1141  tr5 = ti11*cr5+ti12*cr4
1142  ti5 = ti11*ci5+ti12*ci4
1143  tr4 = ti12*cr5-ti11*cr4
1144  ti4 = ti12*ci5-ti11*ci4
1145  ch(i-1,3,k) = tr2+tr5
1146  ch(ic-1,2,k) = tr2-tr5
1147  ch(i,3,k) = ti2+ti5
1148  ch(ic,2,k) = ti5-ti2
1149  ch(i-1,5,k) = tr3+tr4
1150  ch(ic-1,4,k) = tr3-tr4
1151  ch(i,5,k) = ti3+ti4
1152  ch(ic,4,k) = ti4-ti3
1153  102 CONTINUE
1154  103 CONTINUE
1155  RETURN
1156  END
1157 
1158 C> RADFG
1159 C>
1160 C> @param IDO
1161 C> @param IP
1162 C> @param L1
1163 C> @param IDL1
1164 C> @param CC
1165 C> @param C1
1166 C> @param C2
1167 C> @param CH
1168 C> @param CH2
1169 C> @param WA
1170 C>
1171 C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
1172  SUBROUTINE radfg (IDO,IP,L1,IDL1,CC,C1,C2,CH,CH2,WA)
1173  dimension ch(ido,l1,ip) ,cc(ido,ip,l1) ,
1174  1 c1(ido,l1,ip) ,c2(idl1,ip),
1175  2 ch2(idl1,ip) ,wa(1)
1176  DATA tpi/6.28318530717959/
1177  arg = tpi/float(ip)
1178  dcp = cos(arg)
1179  dsp = sin(arg)
1180  ipph = (ip+1)/2
1181  ipp2 = ip+2
1182  idp2 = ido+2
1183  nbd = (ido-1)/2
1184  IF (ido .EQ. 1) GO TO 119
1185  DO 101 ik=1,idl1
1186  ch2(ik,1) = c2(ik,1)
1187  101 CONTINUE
1188  DO 103 j=2,ip
1189  DO 102 k=1,l1
1190  ch(1,k,j) = c1(1,k,j)
1191  102 CONTINUE
1192  103 CONTINUE
1193  IF (nbd .GT. l1) GO TO 107
1194  is = -ido
1195  DO 106 j=2,ip
1196  is = is+ido
1197  idij = is
1198  DO 105 i=3,ido,2
1199  idij = idij+2
1200  DO 104 k=1,l1
1201  ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
1202  ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
1203  104 CONTINUE
1204  105 CONTINUE
1205  106 CONTINUE
1206  GO TO 111
1207  107 is = -ido
1208  DO 110 j=2,ip
1209  is = is+ido
1210  DO 109 k=1,l1
1211  idij = is
1212  DO 108 i=3,ido,2
1213  idij = idij+2
1214  ch(i-1,k,j) = wa(idij-1)*c1(i-1,k,j)+wa(idij)*c1(i,k,j)
1215  ch(i,k,j) = wa(idij-1)*c1(i,k,j)-wa(idij)*c1(i-1,k,j)
1216  108 CONTINUE
1217  109 CONTINUE
1218  110 CONTINUE
1219  111 IF (nbd .LT. l1) GO TO 115
1220  DO 114 j=2,ipph
1221  jc = ipp2-j
1222  DO 113 k=1,l1
1223  DO 112 i=3,ido,2
1224  c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
1225  c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
1226  c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
1227  c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
1228  112 CONTINUE
1229  113 CONTINUE
1230  114 CONTINUE
1231  GO TO 121
1232  115 DO 118 j=2,ipph
1233  jc = ipp2-j
1234  DO 117 i=3,ido,2
1235  DO 116 k=1,l1
1236  c1(i-1,k,j) = ch(i-1,k,j)+ch(i-1,k,jc)
1237  c1(i-1,k,jc) = ch(i,k,j)-ch(i,k,jc)
1238  c1(i,k,j) = ch(i,k,j)+ch(i,k,jc)
1239  c1(i,k,jc) = ch(i-1,k,jc)-ch(i-1,k,j)
1240  116 CONTINUE
1241  117 CONTINUE
1242  118 CONTINUE
1243  GO TO 121
1244  119 DO 120 ik=1,idl1
1245  c2(ik,1) = ch2(ik,1)
1246  120 CONTINUE
1247  121 DO 123 j=2,ipph
1248  jc = ipp2-j
1249  DO 122 k=1,l1
1250  c1(1,k,j) = ch(1,k,j)+ch(1,k,jc)
1251  c1(1,k,jc) = ch(1,k,jc)-ch(1,k,j)
1252  122 CONTINUE
1253  123 CONTINUE
1254 C
1255  ar1 = 1.
1256  ai1 = 0.
1257  DO 127 l=2,ipph
1258  lc = ipp2-l
1259  ar1h = dcp*ar1-dsp*ai1
1260  ai1 = dcp*ai1+dsp*ar1
1261  ar1 = ar1h
1262  DO 124 ik=1,idl1
1263  ch2(ik,l) = c2(ik,1)+ar1*c2(ik,2)
1264  ch2(ik,lc) = ai1*c2(ik,ip)
1265  124 CONTINUE
1266  dc2 = ar1
1267  ds2 = ai1
1268  ar2 = ar1
1269  ai2 = ai1
1270  DO 126 j=3,ipph
1271  jc = ipp2-j
1272  ar2h = dc2*ar2-ds2*ai2
1273  ai2 = dc2*ai2+ds2*ar2
1274  ar2 = ar2h
1275  DO 125 ik=1,idl1
1276  ch2(ik,l) = ch2(ik,l)+ar2*c2(ik,j)
1277  ch2(ik,lc) = ch2(ik,lc)+ai2*c2(ik,jc)
1278  125 CONTINUE
1279  126 CONTINUE
1280  127 CONTINUE
1281  DO 129 j=2,ipph
1282  DO 128 ik=1,idl1
1283  ch2(ik,1) = ch2(ik,1)+c2(ik,j)
1284  128 CONTINUE
1285  129 CONTINUE
1286 C
1287  IF (ido .LT. l1) GO TO 132
1288  DO 131 k=1,l1
1289  DO 130 i=1,ido
1290  cc(i,1,k) = ch(i,k,1)
1291  130 CONTINUE
1292  131 CONTINUE
1293  GO TO 135
1294  132 DO 134 i=1,ido
1295  DO 133 k=1,l1
1296  cc(i,1,k) = ch(i,k,1)
1297  133 CONTINUE
1298  134 CONTINUE
1299  135 DO 137 j=2,ipph
1300  jc = ipp2-j
1301  j2 = j+j
1302  DO 136 k=1,l1
1303  cc(ido,j2-2,k) = ch(1,k,j)
1304  cc(1,j2-1,k) = ch(1,k,jc)
1305  136 CONTINUE
1306  137 CONTINUE
1307  IF (ido .EQ. 1) RETURN
1308  IF (nbd .LT. l1) GO TO 141
1309  DO 140 j=2,ipph
1310  jc = ipp2-j
1311  j2 = j+j
1312  DO 139 k=1,l1
1313  DO 138 i=3,ido,2
1314  ic = idp2-i
1315  cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1316  cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1317  cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1318  cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1319  138 CONTINUE
1320  139 CONTINUE
1321  140 CONTINUE
1322  RETURN
1323  141 DO 144 j=2,ipph
1324  jc = ipp2-j
1325  j2 = j+j
1326  DO 143 i=3,ido,2
1327  ic = idp2-i
1328  DO 142 k=1,l1
1329  cc(i-1,j2-1,k) = ch(i-1,k,j)+ch(i-1,k,jc)
1330  cc(ic-1,j2-2,k) = ch(i-1,k,j)-ch(i-1,k,jc)
1331  cc(i,j2-1,k) = ch(i,k,j)+ch(i,k,jc)
1332  cc(ic,j2-2,k) = ch(i,k,jc)-ch(i,k,j)
1333  142 CONTINUE
1334  143 CONTINUE
1335  144 CONTINUE
1336  RETURN
1337  END
radb5
subroutine radb5(IDO, L1, CC, CH, WA1, WA2, WA3, WA4)
RADB5.
Definition: fftpack.F:691
scfft
subroutine scfft(isign, n, scale, x, y, table, work, isys)
scfft
Definition: fftpack.F:250
radb4
subroutine radb4(IDO, L1, CC, CH, WA1, WA2, WA3)
RADB4.
Definition: fftpack.F:620
rffti
subroutine rffti(N, WSAVE)
RFFTI.
Definition: fftpack.F:310
radf3
subroutine radf3(IDO, L1, CC, CH, WA1, WA2)
RADF3.
Definition: fftpack.F:985
drcft
subroutine drcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
drcft
Definition: fftpack.F:160
radb3
subroutine radb3(IDO, L1, CC, CH, WA1, WA2)
RADB3
Definition: fftpack.F:570
radf5
subroutine radf5(IDO, L1, CC, CH, WA1, WA2, WA3, WA4)
RADF5.
Definition: fftpack.F:1099
radf2
subroutine radf2(IDO, L1, CC, CH, WA1)
RADBG.
Definition: fftpack.F:947
dcrft
subroutine dcrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
dcrft
Definition: fftpack.F:34
rfftb
subroutine rfftb(N, R, WSAVE)
RFFTB.
Definition: fftpack.F:297
scrft
subroutine scrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
scrft
Definition: fftpack.F:80
rfftf1
subroutine rfftf1(N, C, CH, WA, IFAC)
RFFTF1.
Definition: fftpack.F:395
srcft
subroutine srcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
srcft
Definition: fftpack.F:209
csfft
subroutine csfft(isign, n, scale, x, y, table, work, isys)
csfft
Definition: fftpack.F:117
radfg
subroutine radfg(IDO, IP, L1, IDL1, CC, C1, C2, CH, CH2, WA)
RADFG.
Definition: fftpack.F:1173
rffti1
subroutine rffti1(N, WA, IFAC)
RFFTI1.
Definition: fftpack.F:462
radb2
subroutine radb2(IDO, L1, CC, CH, WA1)
RADB2.
Definition: fftpack.F:531
radf4
subroutine radf4(IDO, L1, CC, CH, WA1, WA2, WA3)
RADF4.
Definition: fftpack.F:1032
radbg
subroutine radbg(IDO, IP, L1, IDL1, CC, C1, C2, CH, CH2, WA)
RADBG.
Definition: fftpack.F:769
rfftb1
subroutine rfftb1(N, C, CH, WA, IFAC)
RFFTB1.
Definition: fftpack.F:326
rfftf
subroutine rfftf(N, R, WSAVE)
RFFTF.
Definition: fftpack.F:283