SUBROUTINE FOURIER(A,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN) 8,6FOURIE2B.23
implicit none FOURIE2B.24
! FOURIE2B.25
! Description: FOURIE2B.26
! FOURIE2B.27
! SUBROUTINE 'FOURIER' - MULTIPLE FAST REAL PERIODIC TRANSFORM FOURIE2B.28
! UNIFIED MODEL RE-WRITE OF ECMWF ROUTINE FFT991 FOURIE2B.29
! FOURIE2B.30
! REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT FOURIE2B.31
! OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N FOURIE2B.32
! FOURIE2B.33
! INPUT INFORMATION ... FOURIE2B.34
! A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA FOURIE2B.35
! TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES FOURIE2B.36
! IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N FOURIE2B.37
! INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' FOURIE2B.38
! (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) FOURIE2B.39
! JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR FOURIE2B.40
! N IS THE LENGTH OF THE DATA VECTORS FOURIE2B.41
! LOT IS THE NUMBER OF DATA VECTORS FOURIE2B.42
! ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT FOURIE2B.43
! = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL FOURIE2B.44
! FOURIE2B.45
! ORDERING OF COEFFICIENTS: FOURIE2B.46
! A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) FOURIE2B.47
! WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED FOURIE2B.48
! FOURIE2B.49
! ORDERING OF DATA: FOURIE2B.50
! X(0),X(1),X(2),...,X(N-1), 0 , 0 ; (N+2) LOCATIONS REQUIRED FOURIE2B.51
! FOURIE2B.52
! VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS FOURIE2B.53
! IN PARALLEL FOURIE2B.54
! FOURIE2B.55
! N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE FOURIE2B.56
! EVEN FOURIE2B.57
! FOURIE2B.58
! DEFINITION OF TRANSFORMS: FOURIE2B.59
! ------------------------- FOURIE2B.60
! FOURIE2B.61
! ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) FOURIE2B.62
! WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) FOURIE2B.63
! FOURIE2B.64
! ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) FOURIE2B.65
! B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) FOURIE2B.66
! .................................................................. FOURIE2B.67
! MULTITASKING NOTE: ERROR TRAPPING NOT WORKING IN PRESENT FORM IF FOURIE2B.68
! CODE IS MULTITASKED. FOURIE2B.69
! FOURIE2B.70
! NOT SUITABLE FOR SINGLE COLUMN USE. FOURIE2B.71
! FOURIE2B.72
! VERSION FOR CRAY Y-MP FOURIE2B.73
! REWRITTEN TO UNIFIED MODEL PROGRAMMING STANDARDS FROM ECMWF FOURIE2B.74
! ORIGINAL CODE BY M.H.MAWSON FOURIE2B.75
! ORIGINAL CODE WRITER C. TEMPERTON FOURIE2B.76
! FOURIE2B.77
! History: FOURIE2B.78
! Version Date Comment FOURIE2B.79
! ------- ---- ------- FOURIE2B.80
! FOURIE2B.81
! 4.2* 06/12/96 * Not included in the UM until 4.3. FOURIE2B.82
! This version gives exactly the same results FOURIE2B.83
! on the C90 and T3E as those returned in FOURIE2B.84
! libmet.a on the C90. It deals with inc being FOURIE2B.85
! values other than 1. It does not have the FOURIE2B.86
! following error: FOURIE2B.87
! On the grid to spectral transform it returns FOURIE2B.88
! the the B terms as being positive, whereas FOURIE2B.89
! FFT991 returns the negative for the B terms. FOURIE2B.90
! On the inverse transform a change of sign in FOURIE2B.91
! required again if a change is effected on the FOURIE2B.92
! forward transform. FOURIE2B.93
! Supplied by A Mills, Cray. K Rogers FOURIE2B.94
! 4.3 11/4/97 Corrections made to vn4.2 FOURIE2B.95
! Additional routine InitFourier is not to be FOURIE2B.96
! used by this version of Fourier. It has been FOURIE2B.97
! included for use with a later version FOURIE2B.98
! of Fourier which takes advantage of storing the FOURIE2B.99
! factor and trig tables. FOURIE2B.100
! Supplied by A Mills, Cray. I Edmond FOURIE2B.101
! 4.5 09/7/98 Added direct call to 'INC_1_FOURIER' GBCNF405.68
! if 'INC' is unity. INC_1_FOURIER is GBCNF405.69
! essentially FOURIE2A with some additional GBCNF405.70
! checks. This call avoids any data copies. GBCNF405.71
! Author: Bob Carruthers, Cray Research GBCNF405.72
! FOURIE2B.102
! This routine is compatable with routine in FOURIER FOURIE2B.103
! FOURIE2B.104
! It calls the Cray standard routines scfft and csfft. FOURIE2B.105
! FOURIE2B.106
! It must rearrange the data before and after the call to scfft and csff FOURIE2B.107
! The rearrangement of the data is dependent on whether the transform FOURIE2B.108
! is grid to spectral or spectral to grid. FOURIE2B.109
! FOURIE2B.110
! Be careful to understand all of this before changing anything! FOURIE2B.111
! FOURIE2B.112
! Code Description: FOURIE2B.113
! Language: FORTRAN 77 + common extensions. FOURIE2B.114
! This code is written to UMDP3 v6 programming standards. FOURIE2B.115
! FOURIE2B.116
! System component covered: <appropriate code> FOURIE2B.117
! System Task: <appropriate code> FOURIE2B.118
! FOURIE2B.119
! Declarations: FOURIE2B.120
FOURIE2B.121
! 1.0 Subroutine arguments FOURIE2B.122
! 1.1 Scalar arguments with intent(in): FOURIE2B.123
integer FOURIE2B.124
& inc FOURIE2B.125
&,n FOURIE2B.126
&,lot FOURIE2B.127
&,jump FOURIE2B.128
&,isign FOURIE2B.129
FOURIE2B.130
! 1.2 Array arguments with intent(in): FOURIE2B.131
integer ifax(10) FOURIE2B.132
FOURIE2B.133
real trigs(3*n) FOURIE2B.134
FOURIE2B.135
! 1.3 Array arguments with intent(in/out): FOURIE2B.136
real a(1+jump*(lot-1)+inc*(n+1)) FOURIE2B.137
FOURIE2B.138
! 2.1 Local scalars: FOURIE2B.139
integer FOURIE2B.140
& k FOURIE2B.141
& ,i FOURIE2B.142
& ,l_sign FOURIE2B.143
& ,isy FOURIE2B.144
FOURIE2B.145
real FOURIE2B.146
& scale FOURIE2B.147
& ,local FOURIE2B.148
FOURIE2B.149
! 2.2 Local arrays: FOURIE2B.150
real FOURIE2B.151
& l_work ( 4 + 4 * n) ! this will work for both MPP and PVP FOURIE2B.152
& ,l_table (100 + 4 * n) ! this will work for both MPP and PVP FOURIE2B.153
& ,l_data (n+2) FOURIE2B.154
FOURIE2B.155
! 3.0 Local parameters: FOURIE2B.156
integer forward FOURIE2B.157
parameter(forward=1) ! grid to spectral in scfft FOURIE2B.158
integer backward FOURIE2B.159
parameter(backward=-forward) ! spectral to grid in csff FOURIE2B.160
!- End of header FOURIE2B.161
*IF DEF,T3E GBCNF405.73
c GBCNF405.74
c--check if we can call the T3E Fourier code with no data GBCNF405.75
c copies GBCNF405.76
if(inc.eq.1) then GBCNF405.77
call inc_1_fourier
(a,trigs,ifax,inc,jump,n,lot,isign) GBCNF405.78
return GBCNF405.79
endif GBCNF405.80
*ENDIF GBCNF405.81
FOURIE2B.162
isy = 0 FOURIE2B.163
FOURIE2B.164
call scfft( 0, n, scale, l_data, l_data, l_table, l_work, isy) FOURIE2B.165
if (isign.eq.+1) then FOURIE2B.166
l_sign = backward FOURIE2B.167
scale = 1.0 FOURIE2B.168
do k = 1, lot FOURIE2B.169
do i = 1, n+1, 2 FOURIE2B.170
l_data(i) = a( 1+(k-1)*jump +(i-1)*inc) FOURIE2B.171
l_data(i+1) = -a( 1+(k-1)*jump + i *inc) FOURIE2B.172
enddo FOURIE2B.173
FOURIE2B.174
call csfft(l_sign, n, scale, l_data, l_data, l_table, l_work, isy) FOURIE2B.175
l_data(n+1) = 0.0 FOURIE2B.176
l_data(n+2) = 0.0 FOURIE2B.177
do i = 1, n+2 FOURIE2B.178
a( 1+(k-1)*jump + (i-1) *inc) = l_data(i) FOURIE2B.179
enddo FOURIE2B.180
enddo FOURIE2B.181
else FOURIE2B.182
l_sign = forward FOURIE2B.183
scale = 1.0/n FOURIE2B.184
do k = 1, lot FOURIE2B.185
do i = 1, n FOURIE2B.186
l_data(i) = a( 1+(k-1)*jump + (i-1) *inc) FOURIE2B.187
enddo FOURIE2B.188
l_data(n+1) = 0.0 FOURIE2B.189
l_data(n+2) = 0.0 FOURIE2B.190
FOURIE2B.191
call scfft(l_sign, n, scale, l_data, l_data, l_table, l_work, isy) FOURIE2B.192
do i = 1, n+1, 2 FOURIE2B.193
a( 1+(k-1)*jump +(i-1)*inc) = l_data(i) FOURIE2B.194
a( 1+(k-1)*jump + i *inc) = -l_data(i+1) FOURIE2B.195
enddo FOURIE2B.196
enddo FOURIE2B.197
endif FOURIE2B.198
FOURIE2B.199
end FOURIE2B.200
FOURIE2B.201
*IF DEF,T3E GBCNF405.82
! *****************************COPYRIGHT****************************** GBCNF405.83
! (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved. GBCNF405.84
! GBCNF405.85
! Use, duplication or disclosure of this code is subject to the GBCNF405.86
! restrictions as set forth in the contract. GBCNF405.87
! GBCNF405.88
! Meteorological Office GBCNF405.89
! London Road GBCNF405.90
! BRACKNELL GBCNF405.91
! Berkshire UK GBCNF405.92
! RG12 2SZ GBCNF405.93
! GBCNF405.94
! If no contract has been raised with this copy of the code, the use, GBCNF405.95
! duplication or disclosure of it is strictly prohibited. Permission GBCNF405.96
! to do so must first be obtained in writing from the Head of Numerical GBCNF405.97
! Modelling at the above address. GBCNF405.98
! ******************************COPYRIGHT****************************** GBCNF405.99
!+ GBCNF405.100
! GBCNF405.101
! Subroutine Interface: GBCNF405.102
CLL SUBROUTINE INC_1_FOURIER -------------------------------------- GBCNF405.103
CLL GBCNF405.104
CLL PURPOSE: ...................................................... GBCNF405.105
CLL GBCNF405.106
CLL SUBROUTINE 'INC_1_FOURIER' - MULTIPLE FAST REAL PERIODIC TRANSFORM GBCNF405.107
CLL UNIFIED MODEL RE-WRITE OF ECMWF ROUTINE FFT991, with unit GBCNF405.108
CLL increment data which does therefore need copies on the T3E GBCNF405.109
CLL GBCNF405.110
CLL REAL TRANSFORM OF LENGTH N PERFORMED BY REMOVING REDUNDANT GBCNF405.111
CLL OPERATIONS FROM COMPLEX TRANSFORM OF LENGTH N GBCNF405.112
CLL GBCNF405.113
CLL INPUT INFORMATION ... GBCNF405.114
CLL A IS THE ARRAY CONTAINING INPUT & OUTPUT DATA GBCNF405.115
CLL TRIGS IS A PREVIOUSLY PREPARED LIST OF TRIG FUNCTION VALUES GBCNF405.116
CLL IFAX IS A PREVIOUSLY PREPARED LIST OF FACTORS OF N GBCNF405.117
CLL INC IS THE INCREMENT WITHIN EACH DATA 'VECTOR' GBCNF405.118
CLL (E.G. INC=1 FOR CONSECUTIVELY STORED DATA) GBCNF405.119
CLL GBCNF405.120
CLL INC must be unity! GBCNF405.121
CLL GBCNF405.122
CLL JUMP IS THE INCREMENT BETWEEN THE START OF EACH DATA VECTOR GBCNF405.123
CLL N IS THE LENGTH OF THE DATA VECTORS GBCNF405.124
CLL LOT IS THE NUMBER OF DATA VECTORS GBCNF405.125
CLL ISIGN = +1 FOR TRANSFORM FROM SPECTRAL TO GRIDPOINT GBCNF405.126
CLL = -1 FOR TRANSFORM FROM GRIDPOINT TO SPECTRAL GBCNF405.127
CLL GBCNF405.128
CLL ORDERING OF COEFFICIENTS: GBCNF405.129
CLL A(0),B(0),A(1),B(1),A(2),B(2),...,A(N/2),B(N/2) GBCNF405.130
CLL WHERE B(0)=B(N/2)=0; (N+2) LOCATIONS REQUIRED GBCNF405.131
CLL GBCNF405.132
CLL ORDERING OF DATA: GBCNF405.133
CLL X(0),X(1),X(2),...,X(N-1), 0 , 0 ; (N+2) LOCATIONS REQUIRED GBCNF405.134
CLL GBCNF405.135
CLL VECTORIZATION IS ACHIEVED ON CRAY BY DOING THE TRANSFORMS GBCNF405.136
CLL IN PARALLEL GBCNF405.137
CLL GBCNF405.138
CLL N MUST BE COMPOSED OF FACTORS 2,3 & 5 BUT DOES NOT HAVE TO BE EVEN GBCNF405.139
CLL GBCNF405.140
CLL DEFINITION OF TRANSFORMS: GBCNF405.141
CLL ------------------------- GBCNF405.142
CLL GBCNF405.143
CLL ISIGN=+1: X(J)=SUM(K=0,...,N-1)(C(K)*EXP(2*I*J*K*PI/N)) GBCNF405.144
CLL WHERE C(K)=A(K)+I*B(K) AND C(N-K)=A(K)-I*B(K) GBCNF405.145
CLL GBCNF405.146
CLL ISIGN=-1: A(K)=(1/N)*SUM(J=0,...,N-1)(X(J)*COS(2*J*K*PI/N)) GBCNF405.147
CLL B(K)=-(1/N)*SUM(J=0,...,N-1)(X(J)*SIN(2*J*K*PI/N)) GBCNF405.148
CLL .................................................................. GBCNF405.149
CLL GBCNF405.150
CLL Current Code Owner: Bob Carruthers GBCNF405.151
CLL GBCNF405.152
CLL History: GBCNF405.153
CLL Version Date Comment GBCNF405.154
CLL ======= ==== ======= GBCNF405.155
CLL 4.5 22/06/98 Original simplified code which assumes GBCNF405.156
CLL 'INC' is 1. This avoids lots of copies. GBCNF405.157
CLL Author: Bob Carruthers, Cray Research GBCNF405.158
GBCNF405.159
C*L ARGUMENT LIST GBCNF405.160
GBCNF405.161
subroutine inc_1_fourier(a,trigs,ifax,inc,jump,n,lot,isign) 1,2GBCNF405.162
GBCNF405.163
IMPLICIT NONE GBCNF405.164
GBCNF405.165
INTEGER GBCNF405.166
* INC ! IN INCREMENT BETWEEN ELEMENTS OF DATA VECTOR GBCNF405.167
* ,JUMP ! IN INCREMENT BETWEEN START OF EACH DATA VECTOR GBCNF405.168
* ,N ! IN LENGTH OF DATA VECTOR IN GRID-POINT SPACE GBCNF405.169
* ! WITHOUT EXTRA ZEROES GBCNF405.170
* ,LOT ! IN NUMBER OF DATA VECTORS GBCNF405.171
* ,ISIGN ! IN DETERMINES TYPE OF TRNASFORM GBCNF405.172
* ,IFAX(10) ! IN LIST OF FACTORS OF N GBCNF405.173
GBCNF405.174
REAL GBCNF405.175
* TRIGS(N) ! IN TRIGONOMETRICAL FUNCTIONS GBCNF405.176
GBCNF405.177
REAL GBCNF405.178
* A(JUMP*LOT) ! INOUT DATA GBCNF405.179
GBCNF405.180
C* -------------------------------------------------------------- GBCNF405.181
GBCNF405.182
C*L WORKSPACE. 1 ARRAY IS REQUIRED -------------------------------- GBCNF405.183
GBCNF405.184
*CALL AMAXSIZE GBCNF405.185
GBCNF405.186
integer ii, inc_1_row_length GBCNF405.187
c GBCNF405.188
parameter (inc_1_row_length=2*row_length_max+100) GBCNF405.189
c GBCNF405.190
real inc_1_work((2*n+4)*lot), inc_1_table(inc_1_row_length), GBCNF405.191
2 inc_1_scale GBCNF405.192
c GBCNF405.193
integer inc_1_flag, inc_1_isys, old_n GBCNF405.194
data inc_1_flag/0/, inc_1_isys/0/, old_n/0/ GBCNF405.195
save inc_1_flag, inc_1_isys, old_n GBCNF405.196
c GBCNF405.197
common/inc_1_fourier_fft/inc_1_table GBCNF405.198
GBCNF405.199
C* -------------------------------------------------------------- GBCNF405.200
GBCNF405.201
C*L EXTERNAL ROUTINES -------------------------------------------- GBCNF405.202
GBCNF405.203
C* --------------------------------------------------------------- GBCNF405.204
c GBCNF405.205
c--check the storage space, and compute the Fourier Coefficients GBCNF405.206
if(inc_1_flag.eq.0 .or. n.ne.old_n) then GBCNF405.207
if(inc_1_row_length.lt.(100+2*n)) then GBCNF405.208
write(6,'(/''Table for Cray Research FFT Routines '', GBCNF405.209
2 '' is '',i5,'' Words Long, but needs to be '',i5, GBCNF405.210
3 '' Words Long'')') inc_1_row_length, (100+2*n) GBCNF405.211
call abort() GBCNF405.212
endif GBCNF405.213
inc_1_scale=1.0 GBCNF405.214
call scfft(0, n, inc_1_scale, a, a, GBCNF405.215
2 inc_1_table, inc_1_work, inc_1_isys) GBCNF405.216
inc_1_flag=1 GBCNF405.217
old_n=n GBCNF405.218
endif GBCNF405.219
c GBCNF405.220
c--check that 'INC' is unity GBCNF405.221
if(inc.ne.1) then GBCNF405.222
write(6,'(/''This Version of the Fourier Code '', GBCNF405.223
2 '' only works with INC = 1, not '',i5)') inc GBCNF405.224
call abort() GBCNF405.225
endif GBCNF405.226
c GBCNF405.227
if(isign.eq.-1) then GBCNF405.228
do ii=1,lot GBCNF405.229
inc_1_scale=1./n GBCNF405.230
call scfft(+1, n, inc_1_scale, a((ii-1)*jump+1), GBCNF405.231
2 a((ii-1)*jump+1), inc_1_table, inc_1_work, inc_1_isys) GBCNF405.232
enddo GBCNF405.233
endif GBCNF405.234
c GBCNF405.235
if(isign.eq.+1) then GBCNF405.236
inc_1_scale=1. GBCNF405.237
do ii=1,lot GBCNF405.238
call csfft(-1, n, inc_1_scale, a((ii-1)*jump+1), GBCNF405.239
2 a((ii-1)*jump+1), inc_1_table, inc_1_work, inc_1_isys) GBCNF405.240
enddo GBCNF405.241
endif GBCNF405.242
GBCNF405.243
CL END OF ROUTINE FOURIER GBCNF405.244
GBCNF405.245
RETURN GBCNF405.246
END GBCNF405.247
*ENDIF GBCNF405.248
FOURIE2B.202
!+ FOURIE2B.203
! FOURIE2B.204
! Subroutine Interface: FOURIE2B.205
SUBROUTINE InitFourier(N,IFAX,TRIGS) FOURIE2B.206
implicit none FOURIE2B.207
! FOURIE2B.208
! Description: FOURIE2B.209
! FOURIE2B.210
!C SUBROUTINE 'SET66' - COMPUTES FACTORS OF N & TRIGONOMETRIC FOURIE2B.211
!C FUNCTIONS REQUIRED BY STAGGERED SINE TRANSFORM Var_FFT66 FOURIE2B.212
!C AND STAGGERED COSINE TRANSFORM Var_FFT55 FOURIE2B.213
!C FOURIE2B.214
!C AUTHOR: CLIVE TEMPERTON JANUARY 1987 FOURIE2B.215
!C FOURIE2B.216
! modified Andrew Lorenc Oct 1995: FOURIE2B.217
! the "standard" TRIGS for the FFT routines are put in TRIGS(1:N). FOURIE2B.218
! extra TRIGS for the shifted (co)sine transforms are put in FOURIE2B.219
! TRIGS(N+1:N+N). FOURIE2B.220
! extra TRIGS for the sine transforms are put in TRIGS(2N+1:2N+N). FOURIE2B.221
! FOURIE2B.222
! Parent module: VarMod_Fourier FOURIE2B.223
! FOURIE2B.224
! History: FOURIE2B.225
! Version Date Comment FOURIE2B.226
! ------- ---- ------- FOURIE2B.227
! 4.3 11/4/97 Original code. Alistair Mills FOURIE2B.228
! FOURIE2B.229
! Code Description: FOURIE2B.230
! Language: FORTRAN 77 + common extensions. FOURIE2B.231
! This code is written to UMDP3 v6 programming standards. FOURIE2B.232
! FOURIE2B.233
! System component covered: <appropriate code> FOURIE2B.234
! System Task: <appropriate code> FOURIE2B.235
! FOURIE2B.236
! Declarations: FOURIE2B.237
! FOURIE2B.238
! 1.0 Subroutine arguments FOURIE2B.239
! 1.1 Scalar arguments with intent(in): FOURIE2B.240
Integer FOURIE2B.241
& N FOURIE2B.242
FOURIE2B.243
! 1.2 Array arguments with intent(in): FOURIE2B.244
Integer FOURIE2B.245
& IFAX(10) FOURIE2B.246
FOURIE2B.247
! 1.3 Array arguments with intent(out): FOURIE2B.248
Real FOURIE2B.249
& TRIGS(3*N) FOURIE2B.250
FOURIE2B.251
! 2.1 Local scalars: FOURIE2B.252
Integer FOURIE2B.253
& I FOURIE2B.254
&,K FOURIE2B.255
&,L FOURIE2B.256
&,NU FOURIE2B.257
&,NFAX FOURIE2B.258
&,IFAC FOURIE2B.259
FOURIE2B.260
Real FOURIE2B.261
& DEL FOURIE2B.262
&,NIL FOURIE2B.263
&,NHL FOURIE2B.264
&,ANGLE FOURIE2B.265
FOURIE2B.266
! 2.2 Local arrays: FOURIE2B.267
Integer FOURIE2B.268
& JFAX(10) FOURIE2B.269
&,LFAX(7) FOURIE2B.270
FOURIE2B.271
Character*80 ErrorMessage FOURIE2B.272
FOURIE2B.273
Data LFAX/6,8,5,4,3,2,1/ FOURIE2B.274
FOURIE2B.275
!- End of header FOURIE2B.276
FOURIE2B.277
! TRIGS FOR REAL PERIODIC TRANSFORM FOURIE2B.278
! --------------------------------- FOURIE2B.279
DEL=4.0*ASIN(1.0)/FLOAT(N) FOURIE2B.280
NIL=0 FOURIE2B.281
NHL=(N/2)-1 FOURIE2B.282
!DIR$ IVDEP FOURIE2B.283
! Fujitsu vectorization directive GRB0F405.229
!OCL NOVREC GRB0F405.230
DO 10 K=NIL,NHL FOURIE2B.284
ANGLE=FLOAT(K)*DEL FOURIE2B.285
TRIGS(2*K+1)=COS(ANGLE) FOURIE2B.286
TRIGS(2*K+2)=SIN(ANGLE) FOURIE2B.287
10 CONTINUE FOURIE2B.288
! FOURIE2B.289
! EXTRA TRIGS FOR (CO)SINE TRANSFORM FOURIE2B.290
! ---------------------------------- FOURIE2B.291
FOURIE2B.292
DEL=0.5*DEL FOURIE2B.293
DO K=1,NHL FOURIE2B.294
ANGLE=FLOAT(K)*DEL FOURIE2B.295
TRIGS(2*N+K)=SIN(ANGLE) FOURIE2B.296
END DO FOURIE2B.297
FOURIE2B.298
! EXTRA TRIGS FOR SHIFTED (CO)SINE TRANSFORM FOURIE2B.299
! -------------------------------------- FOURIE2B.300
DEL=0.5*DEL FOURIE2B.301
DO 15 K=1,N FOURIE2B.302
ANGLE=FLOAT(K)*DEL FOURIE2B.303
TRIGS(N+K)=SIN(ANGLE) FOURIE2B.304
15 CONTINUE FOURIE2B.305
FOURIE2B.306
! FOURIE2B.307
! NOW FIND FACTORS OF N FOURIE2B.308
! --------------------- FOURIE2B.309
! FOURIE2B.310
! AL 9/10/95. Removed code following Clive Temperton's advice: FOURIE2B.311
! the "fix for small N" section in SET66 was needed for obscure FOURIE2B.312
! historical reasons in the version which called CAL routines. FOURIE2B.313
! It's not necessary in the all-Fortran version, and removing it FOURIE2B.314
! should make things go a bit faster. FOURIE2B.315
! FOURIE2B.316
! FIND FACTORS OF N (8,6,5,4,3,2; ONLY ONE 8 ALLOWED) FOURIE2B.317
! LOOK FOR SIXES FIRST, STORE FACTORS IN DESCENDING ORDER FOURIE2B.318
NU=N FOURIE2B.319
IFAC=6 FOURIE2B.320
K=0 FOURIE2B.321
L=1 FOURIE2B.322
20 CONTINUE FOURIE2B.323
IF (MOD(NU,IFAC).NE.0) GO TO 30 FOURIE2B.324
K=K+1 FOURIE2B.325
JFAX(K)=IFAC FOURIE2B.326
IF (IFAC.NE.8) GO TO 25 FOURIE2B.327
IF (K.EQ.1) GO TO 25 FOURIE2B.328
JFAX(1)=8 FOURIE2B.329
JFAX(K)=6 FOURIE2B.330
25 CONTINUE FOURIE2B.331
NU=NU/IFAC FOURIE2B.332
IF (NU.EQ.1) GO TO 50 FOURIE2B.333
IF (IFAC.NE.8) GO TO 20 FOURIE2B.334
30 CONTINUE FOURIE2B.335
L=L+1 FOURIE2B.336
IFAC=LFAX(L) FOURIE2B.337
IF (IFAC .GT. 1) GO TO 20 FOURIE2B.338
FOURIE2B.339
! Illegal factors: FOURIE2B.340
write (ErrorMessage, '(a,i4,a)') 'N = ', N, FOURIE2B.341
& ' contains illegal factors.' FOURIE2B.342
FOURIE2B.343
GOTO 9 FOURIE2B.344
! FOURIE2B.345
! NOW REVERSE ORDER OF FACTORS FOURIE2B.346
50 CONTINUE FOURIE2B.347
NFAX=K FOURIE2B.348
IFAX(1)=NFAX FOURIE2B.349
DO 60 I=1,NFAX FOURIE2B.350
IFAX(NFAX+2-I)=JFAX(I) FOURIE2B.351
60 CONTINUE FOURIE2B.352
FOURIE2B.353
9 continue FOURIE2B.354
FOURIE2B.355
END FOURIE2B.356
*ENDIF FOURIE2B.357