SUBROUTINE FOURIER(A,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) 8,6FOURIE2A.22
FOURIE2A.23
IMPLICIT NONE FOURIE2A.24
! Description: FOURIE2A.25
! FOURIE2A.26
! SUBROUTINE 'FOURIER' - MULTIPLE FAST REAL PERIODIC TRANSFORM FOURIE2A.27
! UNIFIED MODEL RE-WRITE OF ECMWF ROUTINE FFT991 FOURIE2A.28
! FOURIE2A.29
! Method: FOURIE2A.30
! REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT FOURIE2A.31
! OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N FOURIE2A.32
! FOURIE2A.33
! INPUT INFORMATION ... FOURIE2A.34
! A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA FOURIE2A.35
! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES FOURIE2A.36
! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N FOURIE2A.37
! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' FOURIE2A.38
! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) FOURIE2A.39
! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR FOURIE2A.40
! N IS THE LENGTH OF THE DATA VECTORS FOURIE2A.41
! LOT IS THE NUMBER OF DATA VECTORS FOURIE2A.42
! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT FOURIE2A.43
! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL FOURIE2A.44
! FOURIE2A.45
! ORDERING OF COEFFICIENTS: FOURIE2A.46
! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) FOURIE2A.47
! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED FOURIE2A.48
! FOURIE2A.49
! ORDERING OF DATA: FOURIE2A.50
! X(0),X(1),X(2),...,X(N-1), 0 , 0 ; (N+2) LOCATIONS REQUIRED FOURIE2A.51
! FOURIE2A.52
! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS FOURIE2A.53
! IN PARALLEL FOURIE2A.54
! FOURIE2A.55
! N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN FOURIE2A.56
! FOURIE2A.57
! DEFINITION OF TRANSFORMS: FOURIE2A.58
! ------------------------- FOURIE2A.59
! FOURIE2A.60
! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) FOURIE2A.61
! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) FOURIE2A.62
! FOURIE2A.63
! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) FOURIE2A.64
! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) FOURIE2A.65
! .................................................................. FOURIE2A.66
! MULTITASKING NOTE: ERROR TRAPPING NOT WORKING IN PRESENT FORM IF FOURIE2A.67
! CODE IS MULTITASKED. FOURIE2A.68
! FOURIE2A.69
! NOT SUITABLE FOR SINGLE COLUMN USE. FOURIE2A.70
! FOURIE2A.71
! Current code owner: M Mawson FOURIE2A.72
! FOURIE2A.73
! History: FOURIE2A.74
! Version Date Comment FOURIE2A.75
! ======= ==== ======= FOURIE2A.76
! 4.1 June 96 Original code at 4.1. Based on modification by FOURIE2A.77
! Cray Research to access Cray specific ffts FOURIE2A.78
! instead of those developed originally by FOURIE2A.79
! ECMWF. R.Rawlins FOURIE2A.80
! FOURIE2A.81
! Code description: FOURIE2A.82
! FORTRAN 77 + common Fortran 90 extensions. FOURIE2A.83
! Written to UM programming standards version 7. FOURIE2A.84
! DOCUMENTATION: NIL. FOURIE2A.85
! FOURIE2A.86
! System component covered: P142 FOURIE2A.87
! FOURIE2A.88
! Subroutine arguments FOURIE2A.89
FOURIE2A.90
FOURIE2A.91
! Scalar arguments with intent(in): FOURIE2A.92
INTEGER FOURIE2A.93
* INC ! IN INCREMENT BETWEEN ELEMENTS OF DATA VECTOR FOURIE2A.94
* ,JUMP ! IN INCREMENT BETWEEN START OF EACH DATA VECTOR FOURIE2A.95
* ,N ! IN LENGTH OF DATA VECTOR IN GRID-POINT SPACE FOURIE2A.96
* ! WITHOUT EXTRA ZEROES FOURIE2A.97
* ,LOT ! IN NUMBER OF DATA VECTORS FOURIE2A.98
* ,ISIGN ! IN DETERMINES TYPE OF TRANSFORM FOURIE2A.99
FOURIE2A.100
FOURIE2A.101
! Array arguments with intent(in): FOURIE2A.102
INTEGER FOURIE2A.103
* IFAX(10) ! IN LIST OF FACTORS OF N FOURIE2A.104
REAL FOURIE2A.105
* TRIGS(N) ! IN TRIGONOMETRICAL FUNCTIONS FOURIE2A.106
FOURIE2A.107
! Array arguments with intent(out): FOURIE2A.108
REAL FOURIE2A.109
* A(JUMP*LOT) ! INOUT DATA FOURIE2A.110
FOURIE2A.111
! Local parameters: FOURIE2A.112
FOURIE2A.113
real cri_work((2*n+4)*lot), cri_table, cri_scale, work FOURIE2A.114
c FOURIE2A.115
integer cri_flag, cri_fft_table, cri_isys, cri_lot FOURIE2A.116
c FOURIE2A.117
parameter (cri_fft_table=1024, cri_lot=128) FOURIE2A.118
c FOURIE2A.119
common/cri_fourier/cri_flag, cri_isys FOURIE2A.120
c FOURIE2A.121
common/cri_fft/cri_table(cri_fft_table) FOURIE2A.122
c FOURIE2A.123
data cri_flag/0/, cri_isys/0/ FOURIE2A.124
C* -------------------------------------------------------------- FOURIE2A.125
FOURIE2A.126
FOURIE2A.127
FOURIE2A.128
! Local scalars FOURIE2A.129
INTEGER FOURIE2A.130
* NFAX ! NUMBER OF FACTORS FOURIE2A.131
*, NX ! N+1 EXCEPT WHERE N IS ODD THEN HOLDS N FOURIE2A.132
*, NBLOX ! NUMBER OF BLOCKS LOT IS SPLIT INTO FOURIE2A.133
*, NVEXREM ! NUMBER OF DATA VECTORS IN BLOCK IF NOT 64 FOURIE2A.134
*, NB ! DO LOOP COUNTER FOURIE2A.135
*, ISTART ! START ADDRESS FOR A BLOCK FOURIE2A.136
*, NVEX ! NUMBER OF ELEMENTS IN VECTOR FOURIE2A.137
*, IA ! USED TO PASS ISTART TO RPASS AND QPASS FOURIE2A.138
*, KI ! VARIABLE USED FOR ADDRESSING FOURIE2A.139
*, IX ! VARIABLE USED FOR ADDRESSING FOURIE2A.140
*, LA ! VARIABLE USED FOR ADDRESSING FOURIE2A.141
*, IGO ! A CONTROL VARIABLE FOURIE2A.142
*, K ! DO LOOP COUNTER FOURIE2A.143
*, IFAC ! HOLDS CURRENT FACTOR FOURIE2A.144
*, IERR ! HOLDS ERROR STATUS FOURIE2A.145
FOURIE2A.146
! Function and subroutine calls: FOURIE2A.147
EXTERNAL FOURIE2A.148
* scfftm,csfftm FOURIE2A.149
FOURIE2A.150
!- End of Header ------------------------------------------------------ FOURIE2A.151
c FOURIE2A.152
if(cri_flag.eq.0) then FOURIE2A.153
if(cri_fft_table.lt.(100+2*n)) then FOURIE2A.154
write(6,'(/''Table for Cray Research FFT Routines '', FOURIE2A.155
2 '' is '',i5,'' Words Long, but needs to '',i5, FOURIE2A.156
3 '' Words Long'')') cri_fft_table, (100+2*n) FOURIE2A.157
call abort
() FOURIE2A.158
endif FOURIE2A.159
cri_scale=1.0 FOURIE2A.160
call scfftm(0, n, lot, cri_scale, a, jump, a, jump/2, FOURIE2A.161
2 cri_table, cri_work, cri_isys) FOURIE2A.162
cri_flag=1 FOURIE2A.163
endif FOURIE2A.164
c FOURIE2A.165
CL FOURIE2A.166
CL ---------------------------------------------------------------- FOURIE2A.167
CL SECTION 1. SET UP INFORMATION FOR SECTIONS 2 AND 3. FOURIE2A.168
CL ---------------------------------------------------------------- FOURIE2A.169
FOURIE2A.170
C SET NUMBER OF FACTORS AND NX FOURIE2A.171
NFAX=IFAX(1) FOURIE2A.172
NX=N+1 FOURIE2A.173
IF (MOD(N,2).EQ.1) NX=N FOURIE2A.174
CL CALCULATE NUMBER OF BLOCKS OF 'cri_lot' DATA VECTORS ARE TO BE SPLIT FOURIE2A.175
NBLOX=1+(LOT-1)/cri_lot FOURIE2A.176
NVEXREM=LOT-(NBLOX-1)*cri_lot FOURIE2A.177
IF (ISIGN.EQ.1) THEN FOURIE2A.178
FOURIE2A.179
CL FOURIE2A.180
CL ---------------------------------------------------------------- FOURIE2A.181
CL SECTION 2. ISIGN=+1, SPECTRAL TO GRIDPOINT TRANSFORM FOURIE2A.182
CL ---------------------------------------------------------------- FOURIE2A.183
FOURIE2A.184
FOURIE2A.185
C LOOP OVER NUMBER OF BLOCKS. FOURIE2A.186
FOURIE2A.187
DO 200 NB=1,NBLOX FOURIE2A.188
IF(NB.EQ.1)THEN FOURIE2A.189
NVEX=LOT-(NBLOX-1)*cri_lot FOURIE2A.190
ISTART=1 FOURIE2A.191
ELSE FOURIE2A.192
NVEX=cri_lot FOURIE2A.193
ISTART=1+NVEXREM*JUMP+(NB-2)*NVEX*JUMP FOURIE2A.194
ENDIF FOURIE2A.195
IA=ISTART FOURIE2A.196
LA=1 FOURIE2A.197
IGO=+1 FOURIE2A.198
FOURIE2A.199
c--call the Cray Research Optimised FFT Routine FOURIE2A.200
cri_scale=1. FOURIE2A.201
call csfftm(-1, n, nvex, cri_scale, a(ia), jump/2, FOURIE2A.202
2 a(ia), jump, cri_table, cri_work, cri_isys) FOURIE2A.203
FOURIE2A.204
C FILL IN ZEROS AT END FOURIE2A.205
C -------------------- FOURIE2A.206
FOURIE2A.207
IX=ISTART+N*INC FOURIE2A.208
DO 220 K=1,NVEX FOURIE2A.209
A(IX)=0.0 FOURIE2A.210
A(IX+INC)=0.0 FOURIE2A.211
IX=IX+JUMP FOURIE2A.212
220 CONTINUE FOURIE2A.213
FOURIE2A.214
200 CONTINUE FOURIE2A.215
ELSE FOURIE2A.216
FOURIE2A.217
CL FOURIE2A.218
CL ---------------------------------------------------------------- FOURIE2A.219
CL SECTION 3. ISIGN=-1, GRIDPOINT TO SPECTRAL TRANSFORM FOURIE2A.220
CL ---------------------------------------------------------------- FOURIE2A.221
FOURIE2A.222
FOURIE2A.223
DO 300 NB=1,NBLOX FOURIE2A.224
IF(NB.EQ.1)THEN FOURIE2A.225
NVEX=LOT-(NBLOX-1)*cri_lot FOURIE2A.226
ISTART=1 FOURIE2A.227
ELSE FOURIE2A.228
NVEX=cri_lot FOURIE2A.229
ISTART=1+NVEXREM*JUMP+(NB-2)*NVEX*JUMP FOURIE2A.230
ENDIF FOURIE2A.231
IA=ISTART FOURIE2A.232
LA=N FOURIE2A.233
IGO=+1 FOURIE2A.234
KI=1 FOURIE2A.235
c--call the Cray Research Optimised FFT Routine FOURIE2A.236
cri_scale=1./n FOURIE2A.237
call scfftm(+1, n, nvex, cri_scale, a(ia), jump, FOURIE2A.238
2 a(ia), jump/2, cri_table, cri_work, cri_isys) FOURIE2A.239
FOURIE2A.240
FOURIE2A.241
300 CONTINUE FOURIE2A.242
END IF FOURIE2A.243
FOURIE2A.244
RETURN FOURIE2A.245
END FOURIE2A.246
FOURIE2A.247
!- End of subroutine code ------------------------------------------ FOURIE2A.248
*ENDIF FOURIE2A.249