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