V3FIT
fft99b.f
1  SUBROUTINE fft99b(work, a, trigs, inc, jump, n, lot)
2  USE stel_kinds
3  IMPLICIT NONE
4 C-----------------------------------------------
5 C D u m m y A r g u m e n t s
6 C-----------------------------------------------
7  INTEGER inc, jump, n, lot
8  REAL(rprec), DIMENSION(*) :: work, a, trigs
9 C-----------------------------------------------
10 C L o c a l P a r a m e t e r s
11 C-----------------------------------------------
12  REAL(rprec), PARAMETER :: one = 1, zero = 0
13 C-----------------------------------------------
14 C L o c a l V a r i a b l e s
15 C-----------------------------------------------
16  INTEGER::nh,nx,ink,ia,ib,ja,jb,l,iabase,ibbase,jabase,jbbase,k
17  REAL(rprec) :: scale, c, s
18 C-----------------------------------------------
19 c
20 c SUBROUTINE fft99b - postprocessing step for fft99, isign=-1
21 c (gridpoint to spectral transform)
22 c
23  nh=n/2
24  nx=n+1
25  ink=inc+inc
26 c
27 c a(0) and a(n/2)
28  scale = one/n
29  ia=1
30  ib=2
31  ja=1
32  jb=n*inc+1
33 cdir$ ivdep
34  DO 10 l=1,lot
35  a(ja)=scale*(work(ia)+work(ib))
36  a(jb)=scale*(work(ia)-work(ib))
37  a(ja+inc) = zero
38  a(jb+inc) = zero
39  ia=ia+nx
40  ib=ib+nx
41  ja=ja+jump
42  jb=jb+jump
43  10 CONTINUE
44 c
45 c remaining wavenumbers
46  scale=0.5_dp*scale
47  iabase=3
48  ibbase=n-1
49  jabase=2*inc+1
50  jbbase=(n-2)*inc+1
51 c
52  DO 30 k=3,nh,2
53  ia=iabase
54  ib=ibbase
55  ja=jabase
56  jb=jbbase
57  c=trigs(n+k)
58  s=trigs(n+k+1)
59 cdir$ ivdep
60  DO 20 l=1,lot
61  a(ja)=scale*((work(ia)+work(ib))
62  * +(c*(work(ia+1)+work(ib+1))+s*(work(ia)-work(ib))))
63  a(jb)=scale*((work(ia)+work(ib))
64  * -(c*(work(ia+1)+work(ib+1))+s*(work(ia)-work(ib))))
65  a(ja+inc)=scale*((c*(work(ia)-work(ib))-s*(work(ia+1)+work(ib+1)))
66  * +(work(ib+1)-work(ia+1)))
67  a(jb+inc)=scale*((c*(work(ia)-work(ib))-s*(work(ia+1)+work(ib+1)))
68  * -(work(ib+1)-work(ia+1)))
69  ia=ia+nx
70  ib=ib+nx
71  ja=ja+jump
72  jb=jb+jump
73  20 CONTINUE
74  iabase=iabase+2
75  ibbase=ibbase-2
76  jabase=jabase+ink
77  jbbase=jbbase-ink
78  30 CONTINUE
79 c
80  IF (iabase.ne.ibbase) GOTO 50
81 c wavenumber n/4 (IF it exists)
82  ia=iabase
83  ja=jabase
84  scale=2*scale
85 cdir$ ivdep
86  DO 40 l=1,lot
87  a(ja)=scale*work(ia)
88  a(ja+inc)=-scale*work(ia+1)
89  ia=ia+nx
90  ja=ja+jump
91  40 CONTINUE
92 c
93  50 CONTINUE
94 
95  END SUBROUTINE fft99b