Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

source: palm/trunk/SOURCE/temperton_fft.f90 @ 1126

Last change on this file since 1126 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list
  • Property svn:keywords set to Id
File size: 51.3 KB
Line 
1 MODULE temperton_fft
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: temperton_fft.f90 392 2009-09-24 10:39:14Z hoffmann $
11!
12! 258 2009-03-13 12:36:03Z heinze
13! Output of messages replaced by message handling routine.
14!
15! Feb. 2007
16! RCS Log replace by Id keyword, revision history cleaned up
17!
18! Revision 1.2  2003/04/16 12:49:25  raasch
19! Abort in case of illegal factors
20!
21! Revision 1.1  2003/03/12 16:41:59  raasch
22! Initial revision
23!
24!
25! Description:
26! ------------
27! Fast Fourier transformation developed by Clive Temperton, ECMWF.
28!------------------------------------------------------------------------------!
29
30    IMPLICIT NONE
31
32    PRIVATE
33
34    PUBLIC set99, fft991cy
35
36
37    INTEGER           ::  nfax(10)   !  array used by *fft991*.
38    REAL, ALLOCATABLE ::  trig(:)    !  array used by *fft991*.
39
40!
41!-- nfft: maximum length of calls to *fft.
42#if defined( __nec )
43    INTEGER, PARAMETER ::  nfft = 256
44#else
45    INTEGER, PARAMETER ::  nfft =  32
46#endif
47
48    INTEGER, PARAMETER ::  nout =   6  ! standard output stream
49
50CONTAINS
51
52  SUBROUTINE fft991cy(a,work,trigs,ifax,inc,jump,n,lot,isign)
53
54    ! Description:
55    !
56    ! Calls fortran-versions of fft's.
57    !
58    ! Method:
59    !
60    ! Subroutine 'fft991cy' - multiple fast real periodic transform
61    ! supersedes previous routine 'fft991cy'.
62    !
63    ! Real transform of length n performed by removing redundant
64    ! operations from complex transform of length n.
65    !
66    ! a       is the array containing input & output data.
67    ! work    is an area of size (n+1)*min(lot,nfft).
68    ! trigs   is a previously prepared list of trig function values.
69    ! ifax    is a previously prepared list of factors of n.
70    ! inc     is the increment within each data 'vector'
71    !         (e.g. inc=1 for consecutively stored data).
72    ! jump    is the increment between the start of each data vector.
73    ! n       is the length of the data vectors.
74    ! lot     is the number of data vectors.
75    ! isign = +1 for transform from spectral to gridpoint
76    !       = -1 for transform from gridpoint to spectral.
77    !
78    ! ordering of coefficients:
79    ! a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
80    ! where b(0)=b(n/2)=0; (n+2) locations required.
81    !
82    ! ordering of data:
83    ! x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
84    !
85    ! Vectorization is achieved on cray by doing the transforms
86    ! in parallel.
87    !
88    ! n must be composed of factors 2,3 & 5 but does not have to be even.
89    !
90    ! definition of transforms:
91    !
92    ! isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n))
93    ! where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
94
95    ! isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n))
96    ! b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n))
97
98    ! calls fortran-versions of fft's  !!!
99    ! dimension a(n),work(n),trigs(n),ifax(1)
100
101
102    IMPLICIT NONE
103
104    !  Scalar arguments
105    INTEGER :: inc, isign, jump, lot, n
106
107    !  Array arguments
108    REAL :: a(*), trigs(*), work(*)
109    INTEGER :: ifax(*)
110
111    !  Local scalars:
112    INTEGER :: i, ia, ibase, ierr, ifac, igo, ii, istart, ix, iz, j, jbase, jj, &
113         &      k, la, nb, nblox, nfax, nvex, nx
114
115    !  Intrinsic functions
116    INTRINSIC MOD
117
118
119    !  Executable statements
120
121    IF (ifax(10)/=n) CALL set99(trigs,ifax,n)
122    nfax = ifax(1)
123    nx = n + 1
124    IF (MOD(n,2)==1) nx = n
125    nblox = 1 + (lot-1)/nfft
126    nvex = lot - (nblox-1)*nfft
127    IF (isign==-1) GO TO 50
128
129    ! isign=+1, spectral to gridpoint transform
130
131    istart = 1
132    DO nb = 1, nblox
133       ia = istart
134       i = istart
135!DIR$ IVDEP
136!CDIR NODEP
137!OCL NOVREC
138       DO j = 1, nvex
139          a(i+inc) = 0.5*a(i)
140          i = i + jump
141       END DO
142       IF (MOD(n,2)==1) GO TO 10
143       i = istart + n*inc
144       DO j = 1, nvex
145          a(i) = 0.5*a(i)
146          i = i + jump
147       END DO
14810     CONTINUE
149       ia = istart + inc
150       la = 1
151       igo = + 1
152
153       DO k = 1, nfax
154          ifac = ifax(k+1)
155          ierr = -1
156          IF (igo==-1) GO TO 20
157          CALL rpassm(a(ia),a(ia+la*inc),work(1),work(ifac*la+1),trigs,inc,1, &
158               &          jump,nx,nvex,n,ifac,la,ierr)
159          GO TO 30
16020        CONTINUE
161          CALL rpassm(work(1),work(la+1),a(ia),a(ia+ifac*la*inc),trigs,1,inc,nx, &
162               &          jump,nvex,n,ifac,la,ierr)
16330        CONTINUE
164          IF (ierr/=0) GO TO 100
165          la = ifac*la
166          igo = -igo
167          ia = istart
168       END DO
169
170       ! If necessary, copy results back to a
171
172       IF (MOD(nfax,2)==0) GO TO 40
173       ibase = 1
174       jbase = ia
175       DO jj = 1, nvex
176          i = ibase
177          j = jbase
178          DO ii = 1, n
179             a(j) = work(i)
180             i = i + 1
181             j = j + inc
182          END DO
183          ibase = ibase + nx
184          jbase = jbase + jump
185       END DO
18640     CONTINUE
187
188       ! Fill in zeros at end
189
190       ix = istart + n*inc
191!DIR$ IVDEP
192!CDIR NODEP
193!OCL NOVREC
194       DO j = 1, nvex
195          a(ix) = 0.0
196          a(ix+inc) = 0.0
197          ix = ix + jump
198       END DO
199
200       istart = istart + nvex*jump
201       nvex = nfft
202    END DO
203    RETURN
204
205    ! isign=-1, gridpoint to spectral transform
206
20750  CONTINUE
208    istart = 1
209    DO nb = 1, nblox
210       ia = istart
211       la = n
212       igo = + 1
213
214       DO k = 1, nfax
215          ifac = ifax(nfax+2-k)
216          la = la/ifac
217          ierr = -1
218          IF (igo==-1) GO TO 60
219          CALL qpassm(a(ia),a(ia+ifac*la*inc),work(1),work(la+1),trigs,inc,1, &
220               &          jump,nx,nvex,n,ifac,la,ierr)
221          GO TO 70
22260        CONTINUE
223          CALL qpassm(work(1),work(ifac*la+1),a(ia),a(ia+la*inc),trigs,1,inc,nx, &
224               &          jump,nvex,n,ifac,la,ierr)
22570        CONTINUE
226          IF (ierr/=0) GO TO 100
227          igo = -igo
228          ia = istart + inc
229       END DO
230
231       ! If necessary, copy results back to a
232
233       IF (MOD(nfax,2)==0) GO TO 80
234       ibase = 1
235       jbase = ia
236       DO jj = 1, nvex
237          i = ibase
238          j = jbase
239          DO ii = 1, n
240             a(j) = work(i)
241             i = i + 1
242             j = j + inc
243          END DO
244          ibase = ibase + nx
245          jbase = jbase + jump
246       END DO
24780     CONTINUE
248
249       ! Shift a(0) & fill in zero imag parts
250
251       ix = istart
252!DIR$ IVDEP
253!CDIR NODEP
254!OCL NOVREC
255       DO j = 1, nvex
256          a(ix) = a(ix+inc)
257          a(ix+inc) = 0.0
258          ix = ix + jump
259       END DO
260       IF (MOD(n,2)==1) GO TO 90
261       iz = istart + (n+1)*inc
262       DO j = 1, nvex
263          a(iz) = 0.0
264          iz = iz + jump
265       END DO
26690     CONTINUE
267
268       istart = istart + nvex*jump
269       nvex = nfft
270    END DO
271    RETURN
272
273    ! Error messages
274
275100 CONTINUE
276
277    SELECT CASE (ierr)
278    CASE (:-1)
279       WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
280    CASE (0)
281       WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
282    CASE (1:)
283       WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
284    END SELECT
285
286    RETURN
287  END SUBROUTINE fft991cy
288
289  SUBROUTINE qpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
290
291    ! Description:
292    !
293    ! Performs one pass through data as part of
294    ! multiple real fft (fourier analysis) routine.
295    !
296    ! Method:
297    !
298    ! a       is first real input vector
299    !         equivalence b(1) with a(ifac*la*inc1+1)
300    ! c       is first real output vector
301    !         equivalence d(1) with c(la*inc2+1)
302    ! trigs   is a precalculated list of sines & cosines
303    ! inc1    is the addressing increment for a
304    ! inc2    is the addressing increment for c
305    ! inc3    is the increment between input vectors a
306    ! inc4    is the increment between output vectors c
307    ! lot     is the number of vectors
308    ! n       is the length of the vectors
309    ! ifac    is the current factor of n
310    !         la = n/(product of factors used so far)
311    ! ierr    is an error indicator:
312    !         0 - pass completed without error
313    !         1 - lot greater than nfft
314    !         2 - ifac not catered for
315    !         3 - ifac only catered for if la=n/ifac
316    !
317
318    IMPLICIT NONE
319
320    !  Scalar arguments
321    INTEGER :: ierr, ifac, inc1, inc2, inc3, inc4, la, lot, n
322
323    !  Array arguments
324    ! REAL :: a(n),b(n),c(n),d(n),trigs(n)
325    REAL :: a(*), b(*), c(*), d(*), trigs(*)
326
327    !  Local scalars:
328    REAL :: a0, a1, a10, a11, a2, a20, a21, a3, a4, a5, a6, b0, b1, b10, b11, &
329         &      b2, b20, b21, b3, b4, b5, b6, c1, c2, c3, c4, c5, qrt5, s1, s2, s3, s4, &
330         &      s5, sin36, sin45, sin60, sin72, z, zqrt5, zsin36, zsin45, zsin60, &
331         &      zsin72
332    INTEGER :: i, ia, ib, ibad, ibase, ic, id, ie, if, ig, igo, ih, iink, ijk, &
333         &      ijump, j, ja, jb, jbase, jc, jd, je, jf, jink, k, kb, kc, kd, ke, kf, &
334         &      kstop, l, m
335
336    !  Intrinsic functions
337    INTRINSIC REAL, SQRT
338
339    !  Data statements
340    DATA sin36/0.587785252292473/, sin72/0.951056516295154/, &
341         &      qrt5/0.559016994374947/, sin60/0.866025403784437/
342
343
344    !  Executable statements
345
346    m = n/ifac
347    iink = la*inc1
348    jink = la*inc2
349    ijump = (ifac-1)*iink
350    kstop = (n-ifac)/(2*ifac)
351
352    ibad = 1
353    IF (lot>nfft) GO TO 180
354    ibase = 0
355    jbase = 0
356    igo = ifac - 1
357    IF (igo==7) igo = 6
358    ibad = 2
359    IF (igo<1 .OR. igo>6) GO TO 180
360    GO TO (10,40,70,100,130,160) igo
361
362    ! Coding for factor 2
363
36410  CONTINUE
365    ia = 1
366    ib = ia + iink
367    ja = 1
368    jb = ja + (2*m-la)*inc2
369
370    IF (la==m) GO TO 30
371
372    DO l = 1, la
373       i = ibase
374       j = jbase
375!DIR$ IVDEP
376!CDIR NODEP
377!OCL NOVREC
378       DO ijk = 1, lot
379          c(ja+j) = a(ia+i) + a(ib+i)
380          c(jb+j) = a(ia+i) - a(ib+i)
381          i = i + inc3
382          j = j + inc4
383       END DO
384       ibase = ibase + inc1
385       jbase = jbase + inc2
386    END DO
387    ja = ja + jink
388    jink = 2*jink
389    jb = jb - jink
390    ibase = ibase + ijump
391    ijump = 2*ijump + iink
392    IF (ja==jb) GO TO 20
393    DO k = la, kstop, la
394       kb = k + k
395       c1 = trigs(kb+1)
396       s1 = trigs(kb+2)
397       jbase = 0
398       DO l = 1, la
399          i = ibase
400          j = jbase
401!DIR$ IVDEP
402!CDIR NODEP
403!OCL NOVREC
404          DO ijk = 1, lot
405             c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
406             c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
407             d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
408             d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
409             i = i + inc3
410             j = j + inc4
411          END DO
412          ibase = ibase + inc1
413          jbase = jbase + inc2
414       END DO
415       ibase = ibase + ijump
416       ja = ja + jink
417       jb = jb - jink
418    END DO
419    IF (ja>jb) GO TO 170
42020  CONTINUE
421    jbase = 0
422    DO l = 1, la
423       i = ibase
424       j = jbase
425!DIR$ IVDEP
426!CDIR NODEP
427!OCL NOVREC
428       DO ijk = 1, lot
429          c(ja+j) = a(ia+i)
430          d(ja+j) = -a(ib+i)
431          i = i + inc3
432          j = j + inc4
433       END DO
434       ibase = ibase + inc1
435       jbase = jbase + inc2
436    END DO
437
438    GO TO 170
43930  CONTINUE
440    z = 1.0/REAL(n)
441    DO l = 1, la
442       i = ibase
443       j = jbase
444!DIR$ IVDEP
445!CDIR NODEP
446!OCL NOVREC
447       DO ijk = 1, lot
448          c(ja+j) = z*(a(ia+i)+a(ib+i))
449          c(jb+j) = z*(a(ia+i)-a(ib+i))
450          i = i + inc3
451          j = j + inc4
452       END DO
453       ibase = ibase + inc1
454       jbase = jbase + inc2
455    END DO
456    GO TO 170
457
458    ! Coding for factor 3
459
46040  CONTINUE
461    ia = 1
462    ib = ia + iink
463    ic = ib + iink
464    ja = 1
465    jb = ja + (2*m-la)*inc2
466    jc = jb
467
468    IF (la==m) GO TO 60
469
470    DO l = 1, la
471       i = ibase
472       j = jbase
473!DIR$ IVDEP
474!CDIR NODEP
475!OCL NOVREC
476       DO ijk = 1, lot
477          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
478          c(jb+j) = a(ia+i) - 0.5*(a(ib+i)+a(ic+i))
479          d(jb+j) = sin60*(a(ic+i)-a(ib+i))
480          i = i + inc3
481          j = j + inc4
482       END DO
483       ibase = ibase + inc1
484       jbase = jbase + inc2
485    END DO
486    ja = ja + jink
487    jink = 2*jink
488    jb = jb + jink
489    jc = jc - jink
490    ibase = ibase + ijump
491    ijump = 2*ijump + iink
492    IF (ja==jc) GO TO 50
493    DO k = la, kstop, la
494       kb = k + k
495       kc = kb + kb
496       c1 = trigs(kb+1)
497       s1 = trigs(kb+2)
498       c2 = trigs(kc+1)
499       s2 = trigs(kc+2)
500       jbase = 0
501       DO l = 1, la
502          i = ibase
503          j = jbase
504!DIR$ IVDEP
505!CDIR NODEP
506!OCL NOVREC
507          DO ijk = 1, lot
508             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
509             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
510             a2 = a(ia+i) - 0.5*a1
511             b2 = b(ia+i) - 0.5*b1
512             a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
513             b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
514             c(ja+j) = a(ia+i) + a1
515             d(ja+j) = b(ia+i) + b1
516             c(jb+j) = a2 + b3
517             d(jb+j) = b2 - a3
518             c(jc+j) = a2 - b3
519             d(jc+j) = -(b2+a3)
520             i = i + inc3
521             j = j + inc4
522          END DO
523          ibase = ibase + inc1
524          jbase = jbase + inc2
525       END DO
526       ibase = ibase + ijump
527       ja = ja + jink
528       jb = jb + jink
529       jc = jc - jink
530    END DO
531    IF (ja>jc) GO TO 170
53250  CONTINUE
533    jbase = 0
534    DO l = 1, la
535       i = ibase
536       j = jbase
537!DIR$ IVDEP
538!CDIR NODEP
539!OCL NOVREC
540       DO ijk = 1, lot
541          c(ja+j) = a(ia+i) + 0.5*(a(ib+i)-a(ic+i))
542          d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
543          c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
544          i = i + inc3
545          j = j + inc4
546       END DO
547       ibase = ibase + inc1
548       jbase = jbase + inc2
549    END DO
550
551    GO TO 170
55260  CONTINUE
553    z = 1.0/REAL(n)
554    zsin60 = z*sin60
555    DO l = 1, la
556       i = ibase
557       j = jbase
558!DIR$ IVDEP
559!CDIR NODEP
560!OCL NOVREC
561       DO ijk = 1, lot
562          c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
563          c(jb+j) = z*(a(ia+i)-0.5*(a(ib+i)+a(ic+i)))
564          d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
565          i = i + inc3
566          j = j + inc4
567       END DO
568       ibase = ibase + inc1
569       jbase = jbase + inc2
570    END DO
571    GO TO 170
572
573    ! Coding for factor 4
574
57570  CONTINUE
576    ia = 1
577    ib = ia + iink
578    ic = ib + iink
579    id = ic + iink
580    ja = 1
581    jb = ja + (2*m-la)*inc2
582    jc = jb + 2*m*inc2
583    jd = jb
584
585    IF (la==m) GO TO 90
586
587    DO l = 1, la
588       i = ibase
589       j = jbase
590!DIR$ IVDEP
591!CDIR NODEP
592!OCL NOVREC
593       DO ijk = 1, lot
594          c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
595          c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
596          c(jb+j) = a(ia+i) - a(ic+i)
597          d(jb+j) = a(id+i) - a(ib+i)
598          i = i + inc3
599          j = j + inc4
600       END DO
601       ibase = ibase + inc1
602       jbase = jbase + inc2
603    END DO
604    ja = ja + jink
605    jink = 2*jink
606    jb = jb + jink
607    jc = jc - jink
608    jd = jd - jink
609    ibase = ibase + ijump
610    ijump = 2*ijump + iink
611    IF (jb==jc) GO TO 80
612    DO k = la, kstop, la
613       kb = k + k
614       kc = kb + kb
615       kd = kc + kb
616       c1 = trigs(kb+1)
617       s1 = trigs(kb+2)
618       c2 = trigs(kc+1)
619       s2 = trigs(kc+2)
620       c3 = trigs(kd+1)
621       s3 = trigs(kd+2)
622       jbase = 0
623       DO l = 1, la
624          i = ibase
625          j = jbase
626!DIR$ IVDEP
627!CDIR NODEP
628!OCL NOVREC
629          DO ijk = 1, lot
630             a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
631             a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
632             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
633             a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
634             b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
635             b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
636             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
637             b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
638             c(ja+j) = a0 + a1
639             c(jc+j) = a0 - a1
640             d(ja+j) = b0 + b1
641             d(jc+j) = b1 - b0
642             c(jb+j) = a2 + b3
643             c(jd+j) = a2 - b3
644             d(jb+j) = b2 - a3
645             d(jd+j) = -(b2+a3)
646             i = i + inc3
647             j = j + inc4
648          END DO
649          ibase = ibase + inc1
650          jbase = jbase + inc2
651       END DO
652       ibase = ibase + ijump
653       ja = ja + jink
654       jb = jb + jink
655       jc = jc - jink
656       jd = jd - jink
657    END DO
658    IF (jb>jc) GO TO 170
65980  CONTINUE
660    sin45 = SQRT(0.5)
661    jbase = 0
662    DO l = 1, la
663       i = ibase
664       j = jbase
665!DIR$ IVDEP
666!CDIR NODEP
667!OCL NOVREC
668       DO ijk = 1, lot
669          c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
670          c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
671          d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
672          d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
673          i = i + inc3
674          j = j + inc4
675       END DO
676       ibase = ibase + inc1
677       jbase = jbase + inc2
678    END DO
679
680    GO TO 170
68190  CONTINUE
682    z = 1.0/REAL(n)
683    DO l = 1, la
684       i = ibase
685       j = jbase
686!DIR$ IVDEP
687!CDIR NODEP
688!OCL NOVREC
689       DO ijk = 1, lot
690          c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
691          c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
692          c(jb+j) = z*(a(ia+i)-a(ic+i))
693          d(jb+j) = z*(a(id+i)-a(ib+i))
694          i = i + inc3
695          j = j + inc4
696       END DO
697       ibase = ibase + inc1
698       jbase = jbase + inc2
699    END DO
700    GO TO 170
701
702    ! Coding for factor 5
703
704100 CONTINUE
705    ia = 1
706    ib = ia + iink
707    ic = ib + iink
708    id = ic + iink
709    ie = id + iink
710    ja = 1
711    jb = ja + (2*m-la)*inc2
712    jc = jb + 2*m*inc2
713    jd = jc
714    je = jb
715
716    IF (la==m) GO TO 120
717
718    DO l = 1, la
719       i = ibase
720       j = jbase
721!DIR$ IVDEP
722!CDIR NODEP
723!OCL NOVREC
724       DO ijk = 1, lot
725          a1 = a(ib+i) + a(ie+i)
726          a3 = a(ib+i) - a(ie+i)
727          a2 = a(ic+i) + a(id+i)
728          a4 = a(ic+i) - a(id+i)
729          a5 = a(ia+i) - 0.25*(a1+a2)
730          a6 = qrt5*(a1-a2)
731          c(ja+j) = a(ia+i) + (a1+a2)
732          c(jb+j) = a5 + a6
733          c(jc+j) = a5 - a6
734          d(jb+j) = -sin72*a3 - sin36*a4
735          d(jc+j) = -sin36*a3 + sin72*a4
736          i = i + inc3
737          j = j + inc4
738       END DO
739       ibase = ibase + inc1
740       jbase = jbase + inc2
741    END DO
742    ja = ja + jink
743    jink = 2*jink
744    jb = jb + jink
745    jc = jc + jink
746    jd = jd - jink
747    je = je - jink
748    ibase = ibase + ijump
749    ijump = 2*ijump + iink
750    IF (jb==jd) GO TO 110
751    DO k = la, kstop, la
752       kb = k + k
753       kc = kb + kb
754       kd = kc + kb
755       ke = kd + kb
756       c1 = trigs(kb+1)
757       s1 = trigs(kb+2)
758       c2 = trigs(kc+1)
759       s2 = trigs(kc+2)
760       c3 = trigs(kd+1)
761       s3 = trigs(kd+2)
762       c4 = trigs(ke+1)
763       s4 = trigs(ke+2)
764       jbase = 0
765       DO l = 1, la
766          i = ibase
767          j = jbase
768!DIR$ IVDEP
769!CDIR NODEP
770!OCL NOVREC
771          DO ijk = 1, lot
772             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
773             a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
774             a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
775             a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
776             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
777             b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
778             b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
779             b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
780             a5 = a(ia+i) - 0.25*(a1+a2)
781             a6 = qrt5*(a1-a2)
782             b5 = b(ia+i) - 0.25*(b1+b2)
783             b6 = qrt5*(b1-b2)
784             a10 = a5 + a6
785             a20 = a5 - a6
786             b10 = b5 + b6
787             b20 = b5 - b6
788             a11 = sin72*b3 + sin36*b4
789             a21 = sin36*b3 - sin72*b4
790             b11 = sin72*a3 + sin36*a4
791             b21 = sin36*a3 - sin72*a4
792             c(ja+j) = a(ia+i) + (a1+a2)
793             c(jb+j) = a10 + a11
794             c(je+j) = a10 - a11
795             c(jc+j) = a20 + a21
796             c(jd+j) = a20 - a21
797             d(ja+j) = b(ia+i) + (b1+b2)
798             d(jb+j) = b10 - b11
799             d(je+j) = -(b10+b11)
800             d(jc+j) = b20 - b21
801             d(jd+j) = -(b20+b21)
802             i = i + inc3
803             j = j + inc4
804          END DO
805          ibase = ibase + inc1
806          jbase = jbase + inc2
807       END DO
808       ibase = ibase + ijump
809       ja = ja + jink
810       jb = jb + jink
811       jc = jc + jink
812       jd = jd - jink
813       je = je - jink
814    END DO
815    IF (jb>jd) GO TO 170
816110 CONTINUE
817    jbase = 0
818    DO l = 1, la
819       i = ibase
820       j = jbase
821!DIR$ IVDEP
822!CDIR NODEP
823!OCL NOVREC
824       DO ijk = 1, lot
825          a1 = a(ib+i) + a(ie+i)
826          a3 = a(ib+i) - a(ie+i)
827          a2 = a(ic+i) + a(id+i)
828          a4 = a(ic+i) - a(id+i)
829          a5 = a(ia+i) + 0.25*(a3-a4)
830          a6 = qrt5*(a3+a4)
831          c(ja+j) = a5 + a6
832          c(jb+j) = a5 - a6
833          c(jc+j) = a(ia+i) - (a3-a4)
834          d(ja+j) = -sin36*a1 - sin72*a2
835          d(jb+j) = -sin72*a1 + sin36*a2
836          i = i + inc3
837          j = j + inc4
838       END DO
839       ibase = ibase + inc1
840       jbase = jbase + inc2
841    END DO
842
843    GO TO 170
844120 CONTINUE
845    z = 1.0/REAL(n)
846    zqrt5 = z*qrt5
847    zsin36 = z*sin36
848    zsin72 = z*sin72
849    DO l = 1, la
850       i = ibase
851       j = jbase
852!DIR$ IVDEP
853!CDIR NODEP
854!OCL NOVREC
855       DO ijk = 1, lot
856          a1 = a(ib+i) + a(ie+i)
857          a3 = a(ib+i) - a(ie+i)
858          a2 = a(ic+i) + a(id+i)
859          a4 = a(ic+i) - a(id+i)
860          a5 = z*(a(ia+i)-0.25*(a1+a2))
861          a6 = zqrt5*(a1-a2)
862          c(ja+j) = z*(a(ia+i)+(a1+a2))
863          c(jb+j) = a5 + a6
864          c(jc+j) = a5 - a6
865          d(jb+j) = -zsin72*a3 - zsin36*a4
866          d(jc+j) = -zsin36*a3 + zsin72*a4
867          i = i + inc3
868          j = j + inc4
869       END DO
870       ibase = ibase + inc1
871       jbase = jbase + inc2
872    END DO
873    GO TO 170
874
875    ! Coding for factor 6
876
877130 CONTINUE
878    ia = 1
879    ib = ia + iink
880    ic = ib + iink
881    id = ic + iink
882    ie = id + iink
883    if = ie + iink
884    ja = 1
885    jb = ja + (2*m-la)*inc2
886    jc = jb + 2*m*inc2
887    jd = jc + 2*m*inc2
888    je = jc
889    jf = jb
890
891    IF (la==m) GO TO 150
892
893    DO l = 1, la
894       i = ibase
895       j = jbase
896!DIR$ IVDEP
897!CDIR NODEP
898!OCL NOVREC
899       DO ijk = 1, lot
900          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
901          c(ja+j) = (a(ia+i)+a(id+i)) + a11
902          c(jc+j) = (a(ia+i)+a(id+i)-0.5*a11)
903          d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
904          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
905          c(jb+j) = (a(ia+i)-a(id+i)) - 0.5*a11
906          d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
907          c(jd+j) = (a(ia+i)-a(id+i)) + a11
908          i = i + inc3
909          j = j + inc4
910       END DO
911       ibase = ibase + inc1
912       jbase = jbase + inc2
913    END DO
914    ja = ja + jink
915    jink = 2*jink
916    jb = jb + jink
917    jc = jc + jink
918    jd = jd - jink
919    je = je - jink
920    jf = jf - jink
921    ibase = ibase + ijump
922    ijump = 2*ijump + iink
923    IF (jc==jd) GO TO 140
924    DO k = la, kstop, la
925       kb = k + k
926       kc = kb + kb
927       kd = kc + kb
928       ke = kd + kb
929       kf = ke + kb
930       c1 = trigs(kb+1)
931       s1 = trigs(kb+2)
932       c2 = trigs(kc+1)
933       s2 = trigs(kc+2)
934       c3 = trigs(kd+1)
935       s3 = trigs(kd+2)
936       c4 = trigs(ke+1)
937       s4 = trigs(ke+2)
938       c5 = trigs(kf+1)
939       s5 = trigs(kf+2)
940       jbase = 0
941       DO l = 1, la
942          i = ibase
943          j = jbase
944!DIR$ IVDEP
945!CDIR NODEP
946!OCL NOVREC
947          DO ijk = 1, lot
948             a1 = c1*a(ib+i) + s1*b(ib+i)
949             b1 = c1*b(ib+i) - s1*a(ib+i)
950             a2 = c2*a(ic+i) + s2*b(ic+i)
951             b2 = c2*b(ic+i) - s2*a(ic+i)
952             a3 = c3*a(id+i) + s3*b(id+i)
953             b3 = c3*b(id+i) - s3*a(id+i)
954             a4 = c4*a(ie+i) + s4*b(ie+i)
955             b4 = c4*b(ie+i) - s4*a(ie+i)
956             a5 = c5*a(if+i) + s5*b(if+i)
957             b5 = c5*b(if+i) - s5*a(if+i)
958             a11 = (a2+a5) + (a1+a4)
959             a20 = (a(ia+i)+a3) - 0.5*a11
960             a21 = sin60*((a2+a5)-(a1+a4))
961             b11 = (b2+b5) + (b1+b4)
962             b20 = (b(ia+i)+b3) - 0.5*b11
963             b21 = sin60*((b2+b5)-(b1+b4))
964             c(ja+j) = (a(ia+i)+a3) + a11
965             d(ja+j) = (b(ia+i)+b3) + b11
966             c(jc+j) = a20 - b21
967             d(jc+j) = a21 + b20
968             c(je+j) = a20 + b21
969             d(je+j) = a21 - b20
970             a11 = (a2-a5) + (a4-a1)
971             a20 = (a(ia+i)-a3) - 0.5*a11
972             a21 = sin60*((a4-a1)-(a2-a5))
973             b11 = (b5-b2) - (b4-b1)
974             b20 = (b3-b(ia+i)) - 0.5*b11
975             b21 = sin60*((b5-b2)+(b4-b1))
976             c(jb+j) = a20 - b21
977             d(jb+j) = a21 - b20
978             c(jd+j) = a11 + (a(ia+i)-a3)
979             d(jd+j) = b11 + (b3-b(ia+i))
980             c(jf+j) = a20 + b21
981             d(jf+j) = a21 + b20
982             i = i + inc3
983             j = j + inc4
984          END DO
985          ibase = ibase + inc1
986          jbase = jbase + inc2
987       END DO
988       ibase = ibase + ijump
989       ja = ja + jink
990       jb = jb + jink
991       jc = jc + jink
992       jd = jd - jink
993       je = je - jink
994       jf = jf - jink
995    END DO
996    IF (jc>jd) GO TO 170
997140 CONTINUE
998    jbase = 0
999    DO l = 1, la
1000       i = ibase
1001       j = jbase
1002!DIR$ IVDEP
1003!CDIR NODEP
1004!OCL NOVREC
1005       DO ijk = 1, lot
1006          c(ja+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
1007          d(ja+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
1008          c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
1009          d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
1010          c(jc+j) = (a(ia+i)+0.5*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
1011          d(jc+j) = -(a(id+i)+0.5*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
1012          i = i + inc3
1013          j = j + inc4
1014       END DO
1015       ibase = ibase + inc1
1016       jbase = jbase + inc2
1017    END DO
1018
1019    GO TO 170
1020150 CONTINUE
1021    z = 1.0/REAL(n)
1022    zsin60 = z*sin60
1023    DO l = 1, la
1024       i = ibase
1025       j = jbase
1026!DIR$ IVDEP
1027!CDIR NODEP
1028!OCL NOVREC
1029       DO ijk = 1, lot
1030          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
1031          c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
1032          c(jc+j) = z*((a(ia+i)+a(id+i))-0.5*a11)
1033          d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
1034          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
1035          c(jb+j) = z*((a(ia+i)-a(id+i))-0.5*a11)
1036          d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1037          c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
1038          i = i + inc3
1039          j = j + inc4
1040       END DO
1041       ibase = ibase + inc1
1042       jbase = jbase + inc2
1043    END DO
1044    GO TO 170
1045
1046    ! Coding for factor 8
1047
1048160 CONTINUE
1049    ibad = 3
1050    IF (la/=m) GO TO 180
1051    ia = 1
1052    ib = ia + iink
1053    ic = ib + iink
1054    id = ic + iink
1055    ie = id + iink
1056    if = ie + iink
1057    ig = if + iink
1058    ih = ig + iink
1059    ja = 1
1060    jb = ja + la*inc2
1061    jc = jb + 2*m*inc2
1062    jd = jc + 2*m*inc2
1063    je = jd + 2*m*inc2
1064    z = 1.0/REAL(n)
1065    zsin45 = z*SQRT(0.5)
1066
1067    DO l = 1, la
1068       i = ibase
1069       j = jbase
1070!DIR$ IVDEP
1071!CDIR NODEP
1072!OCL NOVREC
1073       DO ijk = 1, lot
1074          c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ &
1075               &          a(ih+i))+(a(ib+i)+a(if+i))))
1076          c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ &
1077               &          a(ih+i))+(a(ib+i)+a(if+i))))
1078          c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
1079          d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
1080          c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+ &
1081               &          i)-a(ib+i)))
1082          c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+ &
1083               &          i)-a(ib+i)))
1084          d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + &
1085               &          z*(a(ig+i)-a(ic+i))
1086          d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - &
1087               &          z*(a(ig+i)-a(ic+i))
1088          i = i + inc3
1089          j = j + inc4
1090       END DO
1091       ibase = ibase + inc1
1092       jbase = jbase + inc2
1093    END DO
1094
1095    ! Return
1096
1097170 CONTINUE
1098    ibad = 0
1099180 CONTINUE
1100    ierr = ibad
1101    RETURN
1102  END SUBROUTINE qpassm
1103
1104  SUBROUTINE rpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
1105    ! Dimension a(n),b(n),c(n),d(n),trigs(n)
1106
1107    IMPLICIT NONE
1108
1109    !  Scalar arguments
1110    INTEGER :: ierr, ifac, inc1, inc2, inc3, inc4, la, lot, n
1111
1112    !  Array arguments
1113    REAL :: a(*), b(*), c(*), d(*), trigs(*)
1114
1115    !  Local scalars:
1116    REAL :: c1, c2, c3, c4, c5, qqrt5, qrt5, s1, s2, s3, s4, s5, sin36, sin45, &
1117         &      sin60, sin72, ssin36, ssin45, ssin60, ssin72
1118    INTEGER :: i, ia, ib, ibad, ibase, ic, id, ie, if, igo, iink, ijk, j, ja, &
1119         &      jb, jbase, jc, jd, je, jf, jg, jh, jink, jump, k, kb, kc, kd, ke, kf, &
1120         &      kstop, l, m
1121
1122    !  Local arrays:
1123    REAL :: a10(nfft), a11(nfft), a20(nfft), a21(nfft), b10(nfft), b11(nfft), b20(nfft), &
1124         &      b21(nfft)
1125
1126    !  Intrinsic functions
1127    INTRINSIC SQRT
1128
1129    !  Data statements
1130    DATA sin36/0.587785252292473/, sin72/0.951056516295154/, &
1131         &      qrt5/0.559016994374947/, sin60/0.866025403784437/
1132
1133
1134    !  Executable statements
1135
1136    m = n/ifac
1137    iink = la*inc1
1138    jink = la*inc2
1139    jump = (ifac-1)*jink
1140    kstop = (n-ifac)/(2*ifac)
1141
1142    ibad = 1
1143    IF (lot>nfft) GO TO 180
1144    ibase = 0
1145    jbase = 0
1146    igo = ifac - 1
1147    IF (igo==7) igo = 6
1148    ibad = 2
1149    IF (igo<1 .OR. igo>6) GO TO 180
1150    GO TO (10,40,70,100,130,160) igo
1151
1152    ! Coding for factor 2
1153
115410  CONTINUE
1155    ia = 1
1156    ib = ia + (2*m-la)*inc1
1157    ja = 1
1158    jb = ja + jink
1159
1160    IF (la==m) GO TO 30
1161
1162    DO l = 1, la
1163       i = ibase
1164       j = jbase
1165!DIR$ IVDEP
1166!CDIR NODEP
1167!OCL NOVREC
1168       DO ijk = 1, lot
1169          c(ja+j) = a(ia+i) + a(ib+i)
1170          c(jb+j) = a(ia+i) - a(ib+i)
1171          i = i + inc3
1172          j = j + inc4
1173       END DO
1174       ibase = ibase + inc1
1175       jbase = jbase + inc2
1176    END DO
1177    ia = ia + iink
1178    iink = 2*iink
1179    ib = ib - iink
1180    ibase = 0
1181    jbase = jbase + jump
1182    jump = 2*jump + jink
1183    IF (ia==ib) GO TO 20
1184    DO k = la, kstop, la
1185       kb = k + k
1186       c1 = trigs(kb+1)
1187       s1 = trigs(kb+2)
1188       ibase = 0
1189       DO l = 1, la
1190          i = ibase
1191          j = jbase
1192!DIR$ IVDEP
1193!CDIR NODEP
1194!OCL NOVREC
1195          DO ijk = 1, lot
1196             c(ja+j) = a(ia+i) + a(ib+i)
1197             d(ja+j) = b(ia+i) - b(ib+i)
1198             c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
1199             d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
1200             i = i + inc3
1201             j = j + inc4
1202          END DO
1203          ibase = ibase + inc1
1204          jbase = jbase + inc2
1205       END DO
1206       ia = ia + iink
1207       ib = ib - iink
1208       jbase = jbase + jump
1209    END DO
1210    IF (ia>ib) GO TO 170
121120  CONTINUE
1212    ibase = 0
1213    DO l = 1, la
1214       i = ibase
1215       j = jbase
1216!DIR$ IVDEP
1217!CDIR NODEP
1218!OCL NOVREC
1219       DO ijk = 1, lot
1220          c(ja+j) = a(ia+i)
1221          c(jb+j) = -b(ia+i)
1222          i = i + inc3
1223          j = j + inc4
1224       END DO
1225       ibase = ibase + inc1
1226       jbase = jbase + inc2
1227    END DO
1228
1229    GO TO 170
123030  CONTINUE
1231    DO l = 1, la
1232       i = ibase
1233       j = jbase
1234!DIR$ IVDEP
1235!CDIR NODEP
1236!OCL NOVREC
1237       DO ijk = 1, lot
1238          c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
1239          c(jb+j) = 2.0*(a(ia+i)-a(ib+i))
1240          i = i + inc3
1241          j = j + inc4
1242       END DO
1243       ibase = ibase + inc1
1244       jbase = jbase + inc2
1245    END DO
1246    GO TO 170
1247
1248    ! Coding for factor 3
1249
125040  CONTINUE
1251    ia = 1
1252    ib = ia + (2*m-la)*inc1
1253    ic = ib
1254    ja = 1
1255    jb = ja + jink
1256    jc = jb + jink
1257
1258    IF (la==m) GO TO 60
1259
1260    DO l = 1, la
1261       i = ibase
1262       j = jbase
1263!DIR$ IVDEP
1264!CDIR NODEP
1265!OCL NOVREC
1266       DO ijk = 1, lot
1267          c(ja+j) = a(ia+i) + a(ib+i)
1268          c(jb+j) = (a(ia+i)-0.5*a(ib+i)) - (sin60*(b(ib+i)))
1269          c(jc+j) = (a(ia+i)-0.5*a(ib+i)) + (sin60*(b(ib+i)))
1270          i = i + inc3
1271          j = j + inc4
1272       END DO
1273       ibase = ibase + inc1
1274       jbase = jbase + inc2
1275    END DO
1276    ia = ia + iink
1277    iink = 2*iink
1278    ib = ib + iink
1279    ic = ic - iink
1280    jbase = jbase + jump
1281    jump = 2*jump + jink
1282    IF (ia==ic) GO TO 50
1283    DO k = la, kstop, la
1284       kb = k + k
1285       kc = kb + kb
1286       c1 = trigs(kb+1)
1287       s1 = trigs(kb+2)
1288       c2 = trigs(kc+1)
1289       s2 = trigs(kc+2)
1290       ibase = 0
1291       DO l = 1, la
1292          i = ibase
1293          j = jbase
1294!DIR$ IVDEP
1295!CDIR NODEP
1296!OCL NOVREC
1297          DO ijk = 1, lot
1298             c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1299             d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
1300             c(jb+j) = c1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1301                  &            b(ic+i)))) - s1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
1302                  &            a(ic+i))))
1303             d(jb+j) = s1*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
1304                  &            b(ic+i)))) + c1*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
1305                  &            a(ic+i))))
1306             c(jc+j) = c2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
1307                  &            b(ic+i)))) - s2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1308                  &            a(ic+i))))
1309             d(jc+j) = s2*((a(ia+i)-0.5*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
1310                  &            b(ic+i)))) + c2*((b(ia+i)-0.5*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
1311                  &            a(ic+i))))
1312             i = i + inc3
1313             j = j + inc4
1314          END DO
1315          ibase = ibase + inc1
1316          jbase = jbase + inc2
1317       END DO
1318       ia = ia + iink
1319       ib = ib + iink
1320       ic = ic - iink
1321       jbase = jbase + jump
1322    END DO
1323    IF (ia>ic) GO TO 170
132450  CONTINUE
1325    ibase = 0
1326    DO l = 1, la
1327       i = ibase
1328       j = jbase
1329!DIR$ IVDEP
1330!CDIR NODEP
1331!OCL NOVREC
1332       DO ijk = 1, lot
1333          c(ja+j) = a(ia+i) + a(ib+i)
1334          c(jb+j) = (0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1335          c(jc+j) = -(0.5*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
1336          i = i + inc3
1337          j = j + inc4
1338       END DO
1339       ibase = ibase + inc1
1340       jbase = jbase + inc2
1341    END DO
1342
1343    GO TO 170
134460  CONTINUE
1345    ssin60 = 2.0*sin60
1346    DO l = 1, la
1347       i = ibase
1348       j = jbase
1349!DIR$ IVDEP
1350!CDIR NODEP
1351!OCL NOVREC
1352       DO ijk = 1, lot
1353          c(ja+j) = 2.0*(a(ia+i)+a(ib+i))
1354          c(jb+j) = (2.0*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
1355          c(jc+j) = (2.0*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
1356          i = i + inc3
1357          j = j + inc4
1358       END DO
1359       ibase = ibase + inc1
1360       jbase = jbase + inc2
1361    END DO
1362    GO TO 170
1363
1364    ! Coding for factor 4
1365
136670  CONTINUE
1367    ia = 1
1368    ib = ia + (2*m-la)*inc1
1369    ic = ib + 2*m*inc1
1370    id = ib
1371    ja = 1
1372    jb = ja + jink
1373    jc = jb + jink
1374    jd = jc + jink
1375
1376    IF (la==m) GO TO 90
1377
1378    DO l = 1, la
1379       i = ibase
1380       j = jbase
1381!DIR$ IVDEP
1382!CDIR NODEP
1383!OCL NOVREC
1384       DO ijk = 1, lot
1385          c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
1386          c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
1387          c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
1388          c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
1389          i = i + inc3
1390          j = j + inc4
1391       END DO
1392       ibase = ibase + inc1
1393       jbase = jbase + inc2
1394    END DO
1395    ia = ia + iink
1396    iink = 2*iink
1397    ib = ib + iink
1398    ic = ic - iink
1399    id = id - iink
1400    jbase = jbase + jump
1401    jump = 2*jump + jink
1402    IF (ib==ic) GO TO 80
1403    DO k = la, kstop, la
1404       kb = k + k
1405       kc = kb + kb
1406       kd = kc + kb
1407       c1 = trigs(kb+1)
1408       s1 = trigs(kb+2)
1409       c2 = trigs(kc+1)
1410       s2 = trigs(kc+2)
1411       c3 = trigs(kd+1)
1412       s3 = trigs(kd+2)
1413       ibase = 0
1414       DO l = 1, la
1415          i = ibase
1416          j = jbase
1417!DIR$ IVDEP
1418!CDIR NODEP
1419!OCL NOVREC
1420          DO ijk = 1, lot
1421             c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
1422             d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
1423             c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ &
1424                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
1425             d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ &
1426                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
1427             c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ &
1428                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
1429             d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ &
1430                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
1431             c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ &
1432                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
1433             d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ &
1434                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
1435             i = i + inc3
1436             j = j + inc4
1437          END DO
1438          ibase = ibase + inc1
1439          jbase = jbase + inc2
1440       END DO
1441       ia = ia + iink
1442       ib = ib + iink
1443       ic = ic - iink
1444       id = id - iink
1445       jbase = jbase + jump
1446    END DO
1447    IF (ib>ic) GO TO 170
144880  CONTINUE
1449    ibase = 0
1450    sin45 = SQRT(0.5)
1451    DO l = 1, la
1452       i = ibase
1453       j = jbase
1454!DIR$ IVDEP
1455!CDIR NODEP
1456!OCL NOVREC
1457       DO ijk = 1, lot
1458          c(ja+j) = a(ia+i) + a(ib+i)
1459          c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
1460          c(jc+j) = b(ib+i) - b(ia+i)
1461          c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
1462          i = i + inc3
1463          j = j + inc4
1464       END DO
1465       ibase = ibase + inc1
1466       jbase = jbase + inc2
1467    END DO
1468
1469    GO TO 170
147090  CONTINUE
1471    DO l = 1, la
1472       i = ibase
1473       j = jbase
1474!DIR$ IVDEP
1475!CDIR NODEP
1476!OCL NOVREC
1477       DO ijk = 1, lot
1478          c(ja+j) = 2.0*((a(ia+i)+a(ic+i))+a(ib+i))
1479          c(jb+j) = 2.0*((a(ia+i)-a(ic+i))-b(ib+i))
1480          c(jc+j) = 2.0*((a(ia+i)+a(ic+i))-a(ib+i))
1481          c(jd+j) = 2.0*((a(ia+i)-a(ic+i))+b(ib+i))
1482          i = i + inc3
1483          j = j + inc4
1484       END DO
1485       ibase = ibase + inc1
1486       jbase = jbase + inc2
1487    END DO
1488
1489    ! Coding for factor 5
1490
1491    GO TO 170
1492100 CONTINUE
1493    ia = 1
1494    ib = ia + (2*m-la)*inc1
1495    ic = ib + 2*m*inc1
1496    id = ic
1497    ie = ib
1498    ja = 1
1499    jb = ja + jink
1500    jc = jb + jink
1501    jd = jc + jink
1502    je = jd + jink
1503
1504    IF (la==m) GO TO 120
1505
1506    DO l = 1, la
1507       i = ibase
1508       j = jbase
1509!DIR$ IVDEP
1510!CDIR NODEP
1511!OCL NOVREC
1512       DO ijk = 1, lot
1513          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1514          c(jb+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
1515               &          (sin72*b(ib+i)+sin36*b(ic+i))
1516          c(jc+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
1517               &          (sin36*b(ib+i)-sin72*b(ic+i))
1518          c(jd+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
1519               &          (sin36*b(ib+i)-sin72*b(ic+i))
1520          c(je+j) = ((a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
1521               &          (sin72*b(ib+i)+sin36*b(ic+i))
1522          i = i + inc3
1523          j = j + inc4
1524       END DO
1525       ibase = ibase + inc1
1526       jbase = jbase + inc2
1527    END DO
1528    ia = ia + iink
1529    iink = 2*iink
1530    ib = ib + iink
1531    ic = ic + iink
1532    id = id - iink
1533    ie = ie - iink
1534    jbase = jbase + jump
1535    jump = 2*jump + jink
1536    IF (ib==id) GO TO 110
1537    DO k = la, kstop, la
1538       kb = k + k
1539       kc = kb + kb
1540       kd = kc + kb
1541       ke = kd + kb
1542       c1 = trigs(kb+1)
1543       s1 = trigs(kb+2)
1544       c2 = trigs(kc+1)
1545       s2 = trigs(kc+2)
1546       c3 = trigs(kd+1)
1547       s3 = trigs(kd+2)
1548       c4 = trigs(ke+1)
1549       s4 = trigs(ke+2)
1550       ibase = 0
1551       DO l = 1, la
1552          i = ibase
1553          j = jbase
1554!DIR$ IVDEP
1555!CDIR NODEP
1556!OCL NOVREC
1557          DO ijk = 1, lot
1558
1559             a10(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
1560                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1561             a20(ijk) = (a(ia+i)-0.25*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
1562                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1563             b10(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
1564                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1565             b20(ijk) = (b(ia+i)-0.25*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
1566                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1567             a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
1568             a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
1569             b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
1570             b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
1571
1572             c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
1573             d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
1574             c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
1575             d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
1576             c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
1577             d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
1578             c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
1579             d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
1580             c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
1581             d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
1582
1583             i = i + inc3
1584             j = j + inc4
1585          END DO
1586          ibase = ibase + inc1
1587          jbase = jbase + inc2
1588       END DO
1589       ia = ia + iink
1590       ib = ib + iink
1591       ic = ic + iink
1592       id = id - iink
1593       ie = ie - iink
1594       jbase = jbase + jump
1595    END DO
1596    IF (ib>id) GO TO 170
1597110 CONTINUE
1598    ibase = 0
1599    DO l = 1, la
1600       i = ibase
1601       j = jbase
1602!DIR$ IVDEP
1603!CDIR NODEP
1604!OCL NOVREC
1605       DO ijk = 1, lot
1606          c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
1607          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1608               &          (sin36*b(ia+i)+sin72*b(ib+i))
1609          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1610               &          (sin36*b(ia+i)+sin72*b(ib+i))
1611          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1612               &          (sin72*b(ia+i)-sin36*b(ib+i))
1613          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25*(a(ia+i)+a(ib+i))-a(ic+i))) - &
1614               &          (sin72*b(ia+i)-sin36*b(ib+i))
1615          i = i + inc3
1616          j = j + inc4
1617       END DO
1618       ibase = ibase + inc1
1619       jbase = jbase + inc2
1620    END DO
1621
1622    GO TO 170
1623120 CONTINUE
1624    qqrt5 = 2.0*qrt5
1625    ssin36 = 2.0*sin36
1626    ssin72 = 2.0*sin72
1627    DO l = 1, la
1628       i = ibase
1629       j = jbase
1630!DIR$ IVDEP
1631!CDIR NODEP
1632!OCL NOVREC
1633       DO ijk = 1, lot
1634          c(ja+j) = 2.0*(a(ia+i)+(a(ib+i)+a(ic+i)))
1635          c(jb+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1636               &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
1637          c(jc+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
1638               &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
1639          c(jd+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
1640               &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
1641          c(je+j) = (2.0*(a(ia+i)-0.25*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
1642               &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
1643          i = i + inc3
1644          j = j + inc4
1645       END DO
1646       ibase = ibase + inc1
1647       jbase = jbase + inc2
1648    END DO
1649    GO TO 170
1650
1651    ! Coding for factor 6
1652
1653130 CONTINUE
1654    ia = 1
1655    ib = ia + (2*m-la)*inc1
1656    ic = ib + 2*m*inc1
1657    id = ic + 2*m*inc1
1658    ie = ic
1659    if = ib
1660    ja = 1
1661    jb = ja + jink
1662    jc = jb + jink
1663    jd = jc + jink
1664    je = jd + jink
1665    jf = je + jink
1666
1667    IF (la==m) GO TO 150
1668
1669    DO l = 1, la
1670       i = ibase
1671       j = jbase
1672!DIR$ IVDEP
1673!CDIR NODEP
1674!OCL NOVREC
1675       DO ijk = 1, lot
1676          c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
1677          c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
1678          c(jb+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
1679               &          i)+b(ic+i)))
1680          c(jf+j) = ((a(ia+i)-a(id+i))+0.5*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
1681               &          i)+b(ic+i)))
1682          c(jc+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
1683               &          i)-b(ic+i)))
1684          c(je+j) = ((a(ia+i)+a(id+i))-0.5*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
1685               &          i)-b(ic+i)))
1686          i = i + inc3
1687          j = j + inc4
1688       END DO
1689       ibase = ibase + inc1
1690       jbase = jbase + inc2
1691    END DO
1692    ia = ia + iink
1693    iink = 2*iink
1694    ib = ib + iink
1695    ic = ic + iink
1696    id = id - iink
1697    ie = ie - iink
1698    if = if - iink
1699    jbase = jbase + jump
1700    jump = 2*jump + jink
1701    IF (ic==id) GO TO 140
1702    DO k = la, kstop, la
1703       kb = k + k
1704       kc = kb + kb
1705       kd = kc + kb
1706       ke = kd + kb
1707       kf = ke + kb
1708       c1 = trigs(kb+1)
1709       s1 = trigs(kb+2)
1710       c2 = trigs(kc+1)
1711       s2 = trigs(kc+2)
1712       c3 = trigs(kd+1)
1713       s3 = trigs(kd+2)
1714       c4 = trigs(ke+1)
1715       s4 = trigs(ke+2)
1716       c5 = trigs(kf+1)
1717       s5 = trigs(kf+2)
1718       ibase = 0
1719       DO l = 1, la
1720          i = ibase
1721          j = jbase
1722!DIR$ IVDEP
1723!CDIR NODEP
1724!OCL NOVREC
1725          DO ijk = 1, lot
1726
1727             a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
1728             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5*a11(ijk)
1729             a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
1730             b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
1731             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5*b11(ijk)
1732             b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
1733
1734             c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
1735             d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
1736             c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
1737             d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
1738             c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
1739             d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
1740
1741             a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
1742             b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
1743             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5*a11(ijk)
1744             a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1745             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5*b11(ijk)
1746             b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
1747
1748             c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ &
1749                  &            i))-b11(ijk))
1750             d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ &
1751                  &            i))-b11(ijk))
1752             c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
1753             d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
1754             c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
1755             d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
1756
1757             i = i + inc3
1758             j = j + inc4
1759          END DO
1760          ibase = ibase + inc1
1761          jbase = jbase + inc2
1762       END DO
1763       ia = ia + iink
1764       ib = ib + iink
1765       ic = ic + iink
1766       id = id - iink
1767       ie = ie - iink
1768       if = if - iink
1769       jbase = jbase + jump
1770    END DO
1771    IF (ic>id) GO TO 170
1772140 CONTINUE
1773    ibase = 0
1774    DO l = 1, la
1775       i = ibase
1776       j = jbase
1777!DIR$ IVDEP
1778!CDIR NODEP
1779!OCL NOVREC
1780       DO ijk = 1, lot
1781          c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
1782          c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
1783          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
1784          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5*(b(ia+i)+b(ic+i))+b(ib+i))
1785          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
1786          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5*(a(ia+i)+a(ic+i))-a(ib+i))
1787          i = i + inc3
1788          j = j + inc4
1789       END DO
1790       ibase = ibase + inc1
1791       jbase = jbase + inc2
1792    END DO
1793
1794    GO TO 170
1795150 CONTINUE
1796    ssin60 = 2.0*sin60
1797    DO l = 1, la
1798       i = ibase
1799       j = jbase
1800!DIR$ IVDEP
1801!CDIR NODEP
1802!OCL NOVREC
1803       DO ijk = 1, lot
1804          c(ja+j) = (2.0*(a(ia+i)+a(id+i))) + (2.0*(a(ib+i)+a(ic+i)))
1805          c(jd+j) = (2.0*(a(ia+i)-a(id+i))) - (2.0*(a(ib+i)-a(ic+i)))
1806          c(jb+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
1807               &          i)+b(ic+i)))
1808          c(jf+j) = (2.0*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
1809               &          i)+b(ic+i)))
1810          c(jc+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
1811               &          i)-b(ic+i)))
1812          c(je+j) = (2.0*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
1813               &          i)-b(ic+i)))
1814          i = i + inc3
1815          j = j + inc4
1816       END DO
1817       ibase = ibase + inc1
1818       jbase = jbase + inc2
1819    END DO
1820    GO TO 170
1821
1822    ! Coding for factor 8
1823
1824160 CONTINUE
1825    ibad = 3
1826    IF (la/=m) GO TO 180
1827    ia = 1
1828    ib = ia + la*inc1
1829    ic = ib + 2*la*inc1
1830    id = ic + 2*la*inc1
1831    ie = id + 2*la*inc1
1832    ja = 1
1833    jb = ja + jink
1834    jc = jb + jink
1835    jd = jc + jink
1836    je = jd + jink
1837    jf = je + jink
1838    jg = jf + jink
1839    jh = jg + jink
1840    ssin45 = SQRT(2.0)
1841
1842    DO l = 1, la
1843       i = ibase
1844       j = jbase
1845!DIR$ IVDEP
1846!CDIR NODEP
1847!OCL NOVREC
1848       DO ijk = 1, lot
1849          c(ja+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
1850          c(je+j) = 2.0*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
1851          c(jc+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
1852          c(jg+j) = 2.0*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
1853          c(jb+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
1854               &          i))-(b(ib+i)+b(id+i)))
1855          c(jf+j) = 2.0*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
1856               &          i))-(b(ib+i)+b(id+i)))
1857          c(jd+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
1858               &          i))+(b(ib+i)+b(id+i)))
1859          c(jh+j) = 2.0*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
1860               &          i))+(b(ib+i)+b(id+i)))
1861          i = i + inc3
1862          j = j + inc4
1863       END DO
1864       ibase = ibase + inc1
1865       jbase = jbase + inc2
1866    END DO
1867
1868    ! Return
1869
1870170 CONTINUE
1871    ibad = 0
1872180 CONTINUE
1873    ierr = ibad
1874    RETURN
1875  END SUBROUTINE rpassm
1876
1877  SUBROUTINE set99(trigs,ifax,n)
1878
1879    ! Description:
1880    !
1881    ! Computes factors of n & trigonometric functins required by fft99 & fft991cy
1882    !
1883    ! Method:
1884    !
1885    ! Dimension trigs(n),ifax(1),jfax(10),lfax(7)
1886    !
1887    ! subroutine 'set99' - computes factors of n & trigonometric
1888    ! functins required by fft99 & fft991cy
1889
1890
1891    USE control_parameters
1892    USE pegrid
1893
1894    IMPLICIT NONE
1895
1896    !  Scalar arguments
1897    INTEGER :: n
1898
1899    !  Array arguments
1900    REAL :: trigs(*)
1901    INTEGER :: ifax(*)
1902
1903    !  Local scalars:
1904    REAL :: angle, del
1905    INTEGER :: i, ifac, ixxx, k, l, nfax, nhl, nil, nu
1906
1907    !  Local arrays:
1908    INTEGER :: jfax(10), lfax(7)
1909
1910    !  Intrinsic functions
1911    INTRINSIC ASIN, COS, MOD, REAL, SIN
1912
1913    !  Data statements
1914    DATA lfax/6, 8, 5, 4, 3, 2, 1/
1915
1916
1917    !  Executable statements
1918    ixxx = 1
1919
1920    del = 4.0*ASIN(1.0)/REAL(n)
1921    nil = 0
1922    nhl = (n/2) - 1
1923    DO k = nil, nhl
1924       angle = REAL(k)*del
1925       trigs(2*k+1) = COS(angle)
1926       trigs(2*k+2) = SIN(angle)
1927    END DO
1928
1929    ! Find factors of n (8,6,5,4,3,2; only one 8 allowed)
1930    ! Look for sixes first, store factors in descending order
1931    nu = n
1932    ifac = 6
1933    k = 0
1934    l = 1
193510  CONTINUE
1936    IF (MOD(nu,ifac)/=0) GO TO 30
1937    k = k + 1
1938    jfax(k) = ifac
1939    IF (ifac/=8) GO TO 20
1940    IF (k==1) GO TO 20
1941    jfax(1) = 8
1942    jfax(k) = 6
194320  CONTINUE
1944    nu = nu/ifac
1945    IF (nu==1) GO TO 40
1946    IF (ifac/=8) GO TO 10
194730  CONTINUE
1948    l = l + 1
1949    ifac = lfax(l)
1950    IF (ifac>1) GO TO 10
1951
1952!    WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors'
1953    message_string = 'number of gridpoints along x or/and y ' // &
1954                     'contain illegal  factors' //               &
1955                     '&only factors 8,6,5,4,3,2 are allowed' 
1956    CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
1957
1958    RETURN
1959
1960    ! Now reverse order of factors
196140  CONTINUE
1962    nfax = k
1963    ifax(1) = nfax
1964    DO i = 1, nfax
1965       ifax(nfax+2-i) = jfax(i)
1966    END DO
1967    ifax(10) = n
1968    RETURN
1969  END SUBROUTINE set99
1970
1971 END MODULE temperton_fft
Note: See TracBrowser for help on using the repository browser.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |