1C Copyright 1981-2016 ECMWF.
2C
3C This software is licensed under the terms of the Apache Licence
4C Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5C
6C In applying this licence, ECMWF does not waive the privileges and immunities
7C granted to it by virtue of its status as an intergovernmental organisation
8C nor does it submit to any jurisdiction.
9C
10
11      SUBROUTINE FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
12C
13C---->
14C**** FFT99
15C
16C     PURPOSE
17C     _______
18C
19C     Multiple fast real periodic transform.
20C
21C
22C     INTERFACE
23C     _________
24C
25C     CALL FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
26C
27C
28C     Input parameters
29C     ________________
30C
31C     A      - the array containing input data
32C     WORK   - an area of size (N+1)*MIN(LOT,64)
33C     TRIGS  - a previously prepared list of trig function values
34C     IFAX   - a previously prepared list of factors of N
35C     INC    - the increment within each data 'vector'
36C              (e.g. INC=1 for consecutively stored data)
37C     JUMP   - the increment between the start of each data vector
38C     N      - the length of the data vectors
39C     LOT    - the number of data vectors
40C     ISIGN  -  +1 for transform from spectral to gridpoint
41C            -  -1 for transform from gridpoint to spectral
42C
43C     Output parameters
44C     ________________
45C
46C     A      - the array containing output data
47C
48C
49C     Method
50C     ______
51C
52C     Ordering of coefficients:
53C         A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2)
54C         where B(0) = B(N/2)=0; (N+2) locations required
55C
56C     Ordering of data:
57C         X(N-1),X(0),X(1),X(2),...,X(N-1),X(0)
58C         i.e. explicit cyclic continuity; (N+2) locations required
59C
60C     Vectorization is achieved by doing the transforms in parallel
61C
62C     N must be composed of factors 2,3 & 5 but does not have to be even
63C
64C
65C     Real transform of length N performed by removing redundant
66C     operations from complex transform of length N
67C     Definition of transforms:
68C
69C     ISIGN = +1:
70C         X(J) = SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N))
71C         where C(K) = A(K)+I*B(K) and C(N-K) = A(K)-I*B(K)
72C
73C     ISIGN = -1:
74C         A(K) = (1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N))
75C         B(K) = -(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N))
76C
77C     Externals
78C     _________
79C
80C     RPASSM  - Performs one pass through data as part of multiple real
81C               FFT (fourier analysis) routine.
82C     QPASSM  - Performs one pass through data as part of multiple real
83C               FFT (fourier synthesis) routine.
84C
85C
86C     Reference
87C     _________
88C
89C     None.
90C
91C
92C     Comments
93C     ________
94C
95C     Tidy up of code in older version of same routine.
96C
97C
98C     AUTHOR
99C     ______
100C
101C     J.D.Chambers      *ECMWF*      Nov 1996
102C
103C
104C     MODIFICATIONS
105C     _____________
106C
107C     None.
108C
109C----<
110C     _______________________________________________________
111C
112C*    Section 0. Definition of variables.
113C     _______________________________________________________
114C
115      IMPLICIT NONE
116C
117C     Subroutine arguments
118      REAL A, WORK, TRIGS
119      INTEGER IFAX, INC, JUMP, N, LOT, ISIGN
120      DIMENSION A(N),WORK(N),TRIGS(N),IFAX(10)
121C
122C     Local variables
123      INTEGER NFAX, NX, NBLOX, NVEX, ISTART, NB, IA, I, J
124      INTEGER LA, IGO, K, IFAC, IERR, IBASE, JBASE, JJ, II
125      INTEGER IX, IZ
126C
127C ------------------------------------------------------------------
128C     Section 1.   Initialise
129C ------------------------------------------------------------------
130C
131  100 CONTINUE
132C
133C     Ensure factorization of N has been done.
134      IF(IFAX(10).NE.N) CALL SET99(TRIGS,IFAX,N)
135C
136      NFAX = IFAX(1)
137      IF (MOD(N,2).EQ.1) THEN
138        NX = N
139      ELSE
140        NX = N + 1
141      ENDIF
142C
143C     Calculate number of blocks of 64 vectors and number of vectors
144C     'left over'. This remainder is transformed first.
145C
146      NBLOX = 1 + (LOT-1)/64
147      NVEX  = LOT-(NBLOX-1)*64
148C
149      IF (ISIGN.EQ.1) THEN
150C
151C ------------------------------------------------------------------
152C     Section 2.   Spectral to gridpoint transform.
153C ------------------------------------------------------------------
154C
155  200   CONTINUE
156C
157C       Loop through the blocks of vectors
158        ISTART = 1
159        DO 270 NB = 1,NBLOX
160          IA = ISTART
161          I  = ISTART
162          DO 210 J = 1,NVEX
163            A(I+INC) = 0.5*A(I)
164            I        = I + JUMP
165  210     CONTINUE
166C
167          IF (MOD(N,2).NE.1) THEN
168            I = ISTART + N*INC
169            DO 220 J = 1,NVEX
170              A(I) = 0.5*A(I)
171              I    = I + JUMP
172  220       CONTINUE
173          ENDIF
174C
175          IA  = ISTART + INC
176          LA  = 1
177          IGO = 1
178C
179C         Work through the factors
180          DO 230 K = 1,NFAX
181            IFAC = IFAX(K+1)
182            IERR = -1
183            IF (IGO.NE.-1) THEN
184              CALL RPASSM(A(IA),A(IA+LA*INC),WORK(1),WORK(IFAC*LA+1),
185     X          TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
186            ELSE
187              CALL RPASSM(WORK(1),WORK(LA+1),A(IA),A(IA+IFAC*LA*INC),
188     X          TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
189            ENDIF
190C
191            IF (IERR.NE.0) GO TO 950
192C
193            LA  = IFAC*LA
194            IGO = -IGO
195            IA  = ISTART + INC
196  230     CONTINUE
197C
198C         If necessary, copy results back to A
199C
200          IF (MOD(NFAX,2).NE.0) THEN
201            IBASE = 1
202            JBASE = IA
203            DO 250 JJ = 1,NVEX
204              I = IBASE
205              J = JBASE
206              DO 240 II = 1,N
207                A(J) = WORK(I)
208                I    = I + 1
209                J    = J + INC
210  240         CONTINUE
211              IBASE = IBASE + NX
212              JBASE = JBASE + JUMP
213  250       CONTINUE
214          ENDIF
215C
216C         Fill in cyclic boundary values (ie repeat the data vector
217C         end points at opposite end of the vector)
218C
219          IX = ISTART
220          IZ = ISTART + N*INC
221CDIR$ IVDEP
222          DO 260 J = 1,NVEX
223            A(IX)     = A(IZ)
224            A(IZ+INC) = A(IX+INC)
225            IX        = IX + JUMP
226            IZ        = IZ + JUMP
227  260     CONTINUE
228C
229C         Adjust pointers for next block
230          ISTART = ISTART + NVEX*JUMP
231          NVEX   = 64
232  270   CONTINUE
233C
234      ELSE
235C
236C ------------------------------------------------------------------
237C     Section 3.   Gridpoint to spectral transform.
238C ------------------------------------------------------------------
239C
240  300   CONTINUE
241C
242C       Loop through the blocks of vectors
243        ISTART = 1
244        DO 390 NB  =  1,NBLOX
245          IA  = ISTART + INC
246          LA  = N
247          IGO = 1
248C
249          DO 310 K = 1,NFAX
250            IFAC = IFAX(NFAX+2-K)
251            LA   = LA/IFAC
252            IERR = -1
253            IF (IGO.NE.-1) THEN
254              CALL QPASSM(A(IA),A(IA+IFAC*LA*INC),WORK(1),WORK(LA+1),
255     X          TRIGS,INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR)
256            ELSE
257              CALL QPASSM(WORK(1),WORK(IFAC*LA+1),A(IA),A(IA+LA*INC),
258     X          TRIGS,1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR)
259            ENDIF
260            IF (IERR.NE.0) GO TO 950
261            IGO = -IGO
262            IA  = ISTART + INC
263  310     CONTINUE
264C
265C         If necessary, copy results back to A
266C
267          IF (MOD(NFAX,2).NE.0) THEN
268            IBASE = 1
269            JBASE = IA
270            DO 330 JJ = 1,NVEX
271              I = IBASE
272              J = JBASE
273              DO 320 II = 1,N
274                A(J) = WORK(I)
275                I    = I + 1
276                J    = J + INC
277  320         CONTINUE
278              IBASE = IBASE + NX
279              JBASE = JBASE + JUMP
280  330       CONTINUE
281          ENDIF
282C
283C         Shift A(0) and fill in zero imaginary parts
284C
285          IX = ISTART
286          DO 340 J = 1,NVEX
287            A(IX)     = A(IX+INC)
288            A(IX+INC) = 0.0
289            IX        = IX + JUMP
290  340     CONTINUE
291C
292          IF (MOD(N,2).NE.1) THEN
293            IZ = ISTART + (N+1)*INC
294            DO 350 J = 1,NVEX
295              A(IZ) = 0.0
296              IZ    = IZ + JUMP
297  350       CONTINUE
298C
299          ENDIF
300C
301C         Adjust pointers for next block
302          ISTART = ISTART + NVEX*JUMP
303          NVEX   = 64
304  390   CONTINUE
305C
306      ENDIF
307C
308C ------------------------------------------------------------------
309C     Section 9.   Closedown.
310C ------------------------------------------------------------------
311C
312  900   CONTINUE
313C
314      RETURN
315C
316C     Error messages
317C
318  950 CONTINUE
319      GO TO (960,970,980) IERR
320C
321  960 CONTINUE
322      WRITE(*,*) 'Vector length greater than 64, = ', NVEX
323      GO TO 900
324  970 CONTINUE
325      WRITE(*,*) 'Factor not handled, =', IFAC
326      GO TO 900
327  980 CONTINUE
328      WRITE(*,*) 'Factor only handled if LA*IFAC=N. Factor = ', IFAC
329C
330      END
331