NCEPLIBS-ip 5.2.0
Loading...
Searching...
No Matches
fftpack.F
Go to the documentation of this file.
1C> @file
2C> @brief A concatination of the (FFTPACK)[https://netlib.org/fftpack/] library code.
3C>
4C> FFTPACK is a package of Fortran subprograms for the fast Fourier
5C> transform of periodic and other symmetric sequences. It includes
6C> complex, real, sine, cosine, and quarter-wave transforms.
7C>
8C>Reference:
9C>- P.N. Swarztrauber, Vectorizing the FFTs, in Parallel Computations
10C>(G. Rodrigue, ed.), Academic Press, 1982, pp. 51--83.
11C>
12C> @author Paul N. Swarztrauber, National Center for Atmospheric Research, Boulder, CO
13
14C> dcrft
15C>
16C> @param init
17C> @param x
18C> @param ldx
19C> @param y
20C> @param ldy
21C> @param n
22C> @param m
23C> @param isign
24C> @param scale
25C> @param table
26C> @param n1
27C> @param wrk
28C> @param n2
29C> @param z
30C> @param nz
31C> @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
60C> scrft
61C>
62C> @param init
63C> @param x
64C> @param ldx
65C> @param y
66C> @param ldy
67C> @param n
68C> @param m
69C> @param isign
70C> @param scale
71C> @param table
72C> @param n1
73C> @param wrk
74C> @param n2
75C> @param z
76C> @param nz
77C>
78C> @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
108C> csfft
109C>
110C> @param isign
111C> @param n
112C> @param scale
113C> @param x
114C> @param y
115C> @param table
116C> @param work
117C> @param isys
118C>
119C> @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
143C> drcft
144C>
145C> @param init
146C> @param x
147C> @param ldx
148C> @param y
149C> @param ldy
150C> @param n
151C> @param m
152C> @param isign
153C> @param scale
154C> @param table
155C> @param n1
156C> @param wrk
157C> @param n2
158C> @param z
159C> @param nz
160C>
161C> @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.
186C 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
194C> srcft
195C>
196C> @param init
197C> @param x
198C> @param ldx
199C> @param y
200C> @param ldy
201C> @param n
202C> @param m
203C> @param isign
204C> @param scale
205C> @param table
206C> @param n1
207C> @param wrk
208C> @param n2
209C> @param z
210C> @param nz
211C>
212C> @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.
238C 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
245C> scfft
246C>
247C> @param isign
248C> @param n
249C> @param scale
250C> @param x
251C> @param y
252C> @param table
253C> @param work
254C> @param isys
255C>
256C> @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
283C> RFFTF
284C>
285C> @param N
286C> @param R
287C> @param WSAVE
288C>
289C> @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
297C> RFFTB
298C>
299C> @param N
300C> @param R
301C> @param WSAVE
302C>
303C> @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
311C> RFFTI
312C>
313C> @param N
314C> @param WSAVE
315C>
316C> @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
324C> RFFTB1
325C>
326C> @param N
327C> @param C
328C> @param CH
329C> @param WA
330C> @param IFAC
331C>
332C> @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
393C> RFFTF1
394C>
395C> @param N
396C> @param C
397C> @param CH
398C> @param WA
399C> @param IFAC
400C>
401C> @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
462C> RFFTI1
463C>
464C> @param N
465C> @param WA
466C> @param IFAC
467C>
468C> @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
537C> RADB2
538C>
539C> @param IDO
540C> @param L1
541C> @param CC
542C> @param CH
543C> @param WA1
544C>
545C> @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
581C> RADB3
582C>
583C> @param IDO
584C> @param L1
585C> @param CC
586C> @param CH
587C> @param WA1
588C> @param WA2
589C>
590C> @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
630C> RADB4
631C>
632C> @param IDO
633C> @param L1
634C> @param CC
635C> @param CH
636C> @param WA1
637C> @param WA2
638C> @param WA3
639C>
640C> @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
706C> RADB5
707C>
708C> @param IDO
709C> @param L1
710C> @param CC
711C> @param CH
712C> @param WA1
713C> @param WA2
714C> @param WA3
715C> @param WA4
716C>
717C> @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
782C> RADBG
783C>
784C> @param IDO
785C> @param IP
786C> @param L1
787C> @param IDL1
788C> @param CC
789C> @param C1
790C> @param C2
791C> @param CH
792C> @param CH2
793C> @param WA
794C>
795C> @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
965C> RADBG
966C>
967C> @param IDO
968C> @param L1
969C> @param CC
970C> @param CH
971C> @param WA1
972C>
973C> @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
1008C> RADF3
1009C>
1010C> @param IDO
1011C> @param L1
1012C> @param CC
1013C> @param CH
1014C> @param WA1
1015C> @param WA2
1016C>
1017C> @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
1054C> RADF4
1055C>
1056C> @param IDO
1057C> @param L1
1058C> @param CC
1059C> @param CH
1060C> @param WA1
1061C> @param WA2
1062C> @param WA3
1063C>
1064C> @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
1126C> RADF5
1127C>
1128C> @param IDO
1129C> @param L1
1130C> @param CC
1131C> @param CH
1132C> @param WA1
1133C> @param WA2
1134C> @param WA3
1135C> @param WA4
1136C>
1137C> @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
1198C> RADFG
1199C>
1200C> @param IDO
1201C> @param IP
1202C> @param L1
1203C> @param IDL1
1204C> @param CC
1205C> @param C1
1206C> @param C2
1207C> @param CH
1208C> @param CH2
1209C> @param WA
1210C>
1211C> @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
1294C
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
1326C
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
radf5
subroutine radf5(ido, l1, cc, ch, wa1, wa2, wa3, wa4)
RADF5.
Definition fftpack.F:1139
radbg
subroutine radbg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa)
RADBG.
Definition fftpack.F:797
rffti1
subroutine rffti1(n, wa, ifac)
RFFTI1.
Definition fftpack.F:470
scfft
subroutine scfft(isign, n, scale, x, y, table, work, isys)
scfft
Definition fftpack.F:258
rfftb
subroutine rfftb(n, r, wsave)
RFFTB.
Definition fftpack.F:305
rfftb1
subroutine rfftb1(n, c, ch, wa, ifac)
RFFTB1.
Definition fftpack.F:334
drcft
subroutine drcft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
drcft
Definition fftpack.F:164
rfftf1
subroutine rfftf1(n, c, ch, wa, ifac)
RFFTF1.
Definition fftpack.F:403
dcrft
subroutine dcrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
dcrft
Definition fftpack.F:34
radf4
subroutine radf4(ido, l1, cc, ch, wa1, wa2, wa3)
RADF4.
Definition fftpack.F:1066
rfftf
subroutine rfftf(n, r, wsave)
RFFTF.
Definition fftpack.F:291
radf3
subroutine radf3(ido, l1, cc, ch, wa1, wa2)
RADF3.
Definition fftpack.F:1019
radfg
subroutine radfg(ido, ip, l1, idl1, cc, c1, c2, ch, ch2, wa)
RADFG.
Definition fftpack.F:1213
radb5
subroutine radb5(ido, l1, cc, ch, wa1, wa2, wa3, wa4)
RADB5.
Definition fftpack.F:719
radf2
subroutine radf2(ido, l1, cc, ch, wa1)
RADBG.
Definition fftpack.F:975
radb3
subroutine radb3(ido, l1, cc, ch, wa1, wa2)
RADB3
Definition fftpack.F:592
radb2
subroutine radb2(ido, l1, cc, ch, wa1)
RADB2.
Definition fftpack.F:547
scrft
subroutine scrft(init, x, ldx, y, ldy, n, m, isign, scale, table, n1, wrk, n2, z, nz)
scrft
Definition fftpack.F:82
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
rffti
subroutine rffti(n, wsave)
RFFTI.
Definition fftpack.F:318
radb4
subroutine radb4(ido, l1, cc, ch, wa1, wa2, wa3)
RADB4.
Definition fftpack.F:642