V3FIT
fft99a.f
1  SUBROUTINE fft99a(a, work, 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(*) :: a, work, trigs
9 C-----------------------------------------------
10 C L o c a l V a r i a b l e s
11 C-----------------------------------------------
12  INTEGER::nh,nx,ink,ia,ib,ja,jb,l,iabase,ibbase,jabase,jbbase,k
13  REAL(rprec) :: c, s
14 C-----------------------------------------------
15 c
16 c SUBROUTINE fft99a - preprocessing step for fft99, isign=+1
17 c (spectral to gridpoint transform)
18 c
19  nh=n/2
20  nx=n+1
21  ink=inc+inc
22 c
23 c a(0) and a(n/2)
24  ia=1
25  ib=n*inc+1
26  ja=1
27  jb=2
28 cdir$ ivdep
29  DO 10 l=1,lot
30  work(ja)=a(ia)+a(ib)
31  work(jb)=a(ia)-a(ib)
32  ia=ia+jump
33  ib=ib+jump
34  ja=ja+nx
35  jb=jb+nx
36  10 CONTINUE
37 c
38 c remaining wavenumbers
39  iabase=2*inc+1
40  ibbase=(n-2)*inc+1
41  jabase=3
42  jbbase=n-1
43 c
44  DO 30 k=3,nh,2
45  ia=iabase
46  ib=ibbase
47  ja=jabase
48  jb=jbbase
49  c=trigs(n+k)
50  s=trigs(n+k+1)
51 cdir$ ivdep
52  DO 20 l=1,lot
53  work(ja)=(a(ia)+a(ib))-
54  * (s*(a(ia)-a(ib))+c*(a(ia+inc)+a(ib+inc)))
55  work(jb)=(a(ia)+a(ib))+
56  * (s*(a(ia)-a(ib))+c*(a(ia+inc)+a(ib+inc)))
57  work(ja+1)=(c*(a(ia)-a(ib))-s*(a(ia+inc)+a(ib+inc)))+
58  * (a(ia+inc)-a(ib+inc))
59  work(jb+1)=(c*(a(ia)-a(ib))-s*(a(ia+inc)+a(ib+inc)))-
60  * (a(ia+inc)-a(ib+inc))
61  ia=ia+jump
62  ib=ib+jump
63  ja=ja+nx
64  jb=jb+nx
65  20 CONTINUE
66  iabase=iabase+ink
67  ibbase=ibbase-ink
68  jabase=jabase+2
69  jbbase=jbbase-2
70  30 CONTINUE
71 c
72  IF (iabase.ne.ibbase) GOTO 50
73 c wavenumber n/4 (IF it exists)
74  ia=iabase
75  ja=jabase
76 cdir$ ivdep
77  DO 40 l=1,lot
78  work(ja)=2.0*a(ia)
79  work(ja+1)=-2.0*a(ia+inc)
80  ia=ia+jump
81  ja=ja+nx
82  40 CONTINUE
83 c
84  50 CONTINUE
85 
86  END SUBROUTINE fft99a