[4f24798] | 1 | !c Copyright 1981-2016 ECMWF. |
---|
| 2 | !c |
---|
| 3 | !c This software is licensed under the terms of the Apache Licence |
---|
| 4 | !c Version 2.0 which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. |
---|
| 5 | !c |
---|
| 6 | !c In applying this licence, ECMWF does not waive the privileges and immunities |
---|
| 7 | !c granted to it by virtue of its status as an intergovernmental organisation |
---|
| 8 | !c nor does it submit to any jurisdiction. |
---|
| 9 | !c |
---|
| 10 | |
---|
| 11 | SUBROUTINE FFT99(A,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) |
---|
| 12 | DIMENSION A(N),WORK(N),TRIGS(N),IFAX(*) |
---|
| 13 | !C |
---|
| 14 | !C SUBROUTINE 'FFT99' - MULTIPLE FAST REAL PERIODIC TRANSFORM |
---|
| 15 | !C SUPERSEDES PREVIOUS ROUTINE 'FFT99' |
---|
| 16 | !C |
---|
| 17 | !C REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT |
---|
| 18 | !C OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N |
---|
| 19 | !C |
---|
| 20 | !C A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA |
---|
| 21 | !C WORK IS AN AREA OF SIZE (N+1)*MIN(LOT,64) |
---|
| 22 | !C TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES |
---|
| 23 | !C IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N |
---|
| 24 | !C INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' |
---|
| 25 | !C (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) |
---|
| 26 | !C JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR |
---|
| 27 | !C N IS THE LENGTH OF THE DATA VECTORS |
---|
| 28 | !C LOT IS THE NUMBER OF DATA VECTORS |
---|
| 29 | !C ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT |
---|
| 30 | !C = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL |
---|
| 31 | !C |
---|
| 32 | !C ORDERING OF COEFFICIENTS: |
---|
| 33 | !C A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) |
---|
| 34 | !C WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED |
---|
| 35 | !C |
---|
| 36 | !C ORDERING OF DATA: |
---|
| 37 | !C X(N-1),X(0),X(1),X(2),...,X(N-1),X(0) |
---|
| 38 | !C I.E. EXPLICIT CYCLIC CONTINUITY; (N+2) LOCATIONS REQUIRED |
---|
| 39 | !C |
---|
| 40 | !C VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS |
---|
| 41 | !C IN PARALLEL |
---|
| 42 | !C |
---|
| 43 | !C N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN |
---|
| 44 | !C |
---|
| 45 | !C DEFINITION OF TRANSFORMS: |
---|
| 46 | !C ------------------------- |
---|
| 47 | !C |
---|
| 48 | !C ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) |
---|
| 49 | !C WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) |
---|
| 50 | !C |
---|
| 51 | !C ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) |
---|
| 52 | !C B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) |
---|
| 53 | !C |
---|
| 54 | IF(IFAX(10).NE.N) CALL SET99(TRIGS,IFAX,N) |
---|
| 55 | NFAX=IFAX(1) |
---|
| 56 | NX=N+1 |
---|
| 57 | IF (MOD(N,2).EQ.1) NX=N |
---|
| 58 | NBLOX=1+(LOT-1)/64 |
---|
| 59 | NVEX=LOT-(NBLOX-1)*64 |
---|
| 60 | IF (ISIGN.EQ.-1) GO TO 300 |
---|
| 61 | !C |
---|
| 62 | !C ISIGN=+1, SPECTRAL TO GRIDPOINT TRANSFORM |
---|
| 63 | !C ----------------------------------------- |
---|
| 64 | 100 CONTINUE |
---|
| 65 | ISTART=1 |
---|
| 66 | DO 220 NB=1,NBLOX |
---|
| 67 | IA=ISTART |
---|
| 68 | I=ISTART |
---|
| 69 | DO 110 J=1,NVEX |
---|
| 70 | A(I+INC)=0.5*A(I) |
---|
| 71 | I=I+JUMP |
---|
| 72 | 110 CONTINUE |
---|
| 73 | IF (MOD(N,2).EQ.1) GO TO 130 |
---|
| 74 | I=ISTART+N*INC |
---|
| 75 | DO 120 J=1,NVEX |
---|
| 76 | A(I)=0.5*A(I) |
---|
| 77 | I=I+JUMP |
---|
| 78 | 120 CONTINUE |
---|
| 79 | 130 CONTINUE |
---|
| 80 | IA=ISTART+INC |
---|
| 81 | LA=1 |
---|
| 82 | IGO=+1 |
---|
| 83 | !C |
---|
| 84 | DO 160 K=1,NFAX |
---|
| 85 | IFAC=IFAX(K+1) |
---|
| 86 | IERR=-1 |
---|
| 87 | IF (IGO.EQ.-1) GO TO 140 |
---|
| 88 | CALL RPASSM(A(IA),A(IA+LA*INC),WORK(1),WORK(IFAC*LA+1),TRIGS,& |
---|
| 89 | INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR) |
---|
| 90 | GO TO 150 |
---|
| 91 | 140 CONTINUE |
---|
| 92 | CALL RPASSM(WORK(1),WORK(LA+1),A(IA),A(IA+IFAC*LA*INC),TRIGS,& |
---|
| 93 | 1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) |
---|
| 94 | 150 CONTINUE |
---|
| 95 | IF (IERR.NE.0) GO TO 500 |
---|
| 96 | LA=IFAC*LA |
---|
| 97 | IGO=-IGO |
---|
| 98 | IA=ISTART+INC |
---|
| 99 | 160 CONTINUE |
---|
| 100 | !C |
---|
| 101 | !C IF NECESSARY, COPY RESULTS BACK TO A |
---|
| 102 | !C ------------------------------------ |
---|
| 103 | IF (MOD(NFAX,2).EQ.0) GO TO 190 |
---|
| 104 | IBASE=1 |
---|
| 105 | JBASE=IA |
---|
| 106 | DO 180 JJ=1,NVEX |
---|
| 107 | I=IBASE |
---|
| 108 | J=JBASE |
---|
| 109 | DO 170 II=1,N |
---|
| 110 | A(J)=WORK(I) |
---|
| 111 | I=I+1 |
---|
| 112 | J=J+INC |
---|
| 113 | 170 CONTINUE |
---|
| 114 | IBASE=IBASE+NX |
---|
| 115 | JBASE=JBASE+JUMP |
---|
| 116 | 180 CONTINUE |
---|
| 117 | 190 CONTINUE |
---|
| 118 | !C |
---|
| 119 | !C FILL IN CYCLIC BOUNDARY VALUES |
---|
| 120 | !C ------------------------------ |
---|
| 121 | IX=ISTART |
---|
| 122 | IZ=ISTART+N*INC |
---|
| 123 | !DIR$ IVDEP |
---|
| 124 | DO 200 J=1,NVEX |
---|
| 125 | A(IX)=A(IZ) |
---|
| 126 | A(IZ+INC)=A(IX+INC) |
---|
| 127 | IX=IX+JUMP |
---|
| 128 | IZ=IZ+JUMP |
---|
| 129 | 200 CONTINUE |
---|
| 130 | !C |
---|
| 131 | ISTART=ISTART+NVEX*JUMP |
---|
| 132 | NVEX=64 |
---|
| 133 | 220 CONTINUE |
---|
| 134 | RETURN |
---|
| 135 | !C |
---|
| 136 | !C ISIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM |
---|
| 137 | !C ----------------------------------------- |
---|
| 138 | 300 CONTINUE |
---|
| 139 | ISTART=1 |
---|
| 140 | DO 410 NB=1,NBLOX |
---|
| 141 | IA=ISTART+INC |
---|
| 142 | LA=N |
---|
| 143 | IGO=+1 |
---|
| 144 | !C |
---|
| 145 | DO 340 K=1,NFAX |
---|
| 146 | IFAC=IFAX(NFAX+2-K) |
---|
| 147 | LA=LA/IFAC |
---|
| 148 | IERR=-1 |
---|
| 149 | IF (IGO.EQ.-1) GO TO 320 |
---|
| 150 | CALL QPASSM(A(IA),A(IA+IFAC*LA*INC),WORK(1),WORK(LA+1),TRIGS,& |
---|
| 151 | INC,1,JUMP,NX,NVEX,N,IFAC,LA,IERR) |
---|
| 152 | GO TO 330 |
---|
| 153 | 320 CONTINUE |
---|
| 154 | CALL QPASSM(WORK(1),WORK(IFAC*LA+1),A(IA),A(IA+LA*INC),TRIGS,& |
---|
| 155 | 1,INC,NX,JUMP,NVEX,N,IFAC,LA,IERR) |
---|
| 156 | 330 CONTINUE |
---|
| 157 | IF (IERR.NE.0) GO TO 500 |
---|
| 158 | IGO=-IGO |
---|
| 159 | IA=ISTART+INC |
---|
| 160 | 340 CONTINUE |
---|
| 161 | !C |
---|
| 162 | !C IF NECESSARY, COPY RESULTS BACK TO A |
---|
| 163 | !C ------------------------------------ |
---|
| 164 | IF (MOD(NFAX,2).EQ.0) GO TO 370 |
---|
| 165 | IBASE=1 |
---|
| 166 | JBASE=IA |
---|
| 167 | DO 360 JJ=1,NVEX |
---|
| 168 | I=IBASE |
---|
| 169 | J=JBASE |
---|
| 170 | DO 350 II=1,N |
---|
| 171 | A(J)=WORK(I) |
---|
| 172 | I=I+1 |
---|
| 173 | J=J+INC |
---|
| 174 | 350 CONTINUE |
---|
| 175 | IBASE=IBASE+NX |
---|
| 176 | JBASE=JBASE+JUMP |
---|
| 177 | 360 CONTINUE |
---|
| 178 | 370 CONTINUE |
---|
| 179 | !C |
---|
| 180 | !C SHIFT A(0) & FILL IN ZERO IMAG PARTS |
---|
| 181 | !C ------------------------------------ |
---|
| 182 | IX=ISTART |
---|
| 183 | DO 380 J=1,NVEX |
---|
| 184 | A(IX)=A(IX+INC) |
---|
| 185 | A(IX+INC)=0.0 |
---|
| 186 | IX=IX+JUMP |
---|
| 187 | 380 CONTINUE |
---|
| 188 | IF (MOD(N,2).EQ.1) GO TO 400 |
---|
| 189 | IZ=ISTART+(N+1)*INC |
---|
| 190 | DO 390 J=1,NVEX |
---|
| 191 | A(IZ)=0.0 |
---|
| 192 | IZ=IZ+JUMP |
---|
| 193 | 390 CONTINUE |
---|
| 194 | 400 CONTINUE |
---|
| 195 | !C |
---|
| 196 | ISTART=ISTART+NVEX*JUMP |
---|
| 197 | NVEX=64 |
---|
| 198 | 410 CONTINUE |
---|
| 199 | RETURN |
---|
| 200 | !C |
---|
| 201 | !C ERROR MESSAGES |
---|
| 202 | !C -------------- |
---|
| 203 | 500 CONTINUE |
---|
| 204 | GO TO (510,530,550) IERR |
---|
| 205 | 510 CONTINUE |
---|
| 206 | WRITE(6,520) NVEX |
---|
| 207 | 520 FORMAT('1VECTOR LENGTH =',I4,', GREATER THAN 64') |
---|
| 208 | GO TO 570 |
---|
| 209 | 530 CONTINUE |
---|
| 210 | WRITE(6,540) IFAC |
---|
| 211 | 540 FORMAT( '1FACTOR =',I3,', NOT CATERED FOR') |
---|
| 212 | GO TO 570 |
---|
| 213 | 550 CONTINUE |
---|
| 214 | WRITE(6,560) IFAC |
---|
| 215 | 560 FORMAT('1FACTOR =',I3,', ONLY CATERED FOR IF LA*IFAC=N') |
---|
| 216 | 570 CONTINUE |
---|
| 217 | RETURN |
---|
| 218 | END |
---|