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