source: palm/trunk/SOURCE/temperton_fft_mod.f90

Last change on this file was 4828, checked in by Giersch, 4 months ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

  • Property svn:keywords set to Id
File size: 141.8 KB
Line 
1!> @file temperton_fft_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2021 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: temperton_fft_mod.f90 4828 2021-01-05 11:21:41Z banzhafs $
27! File re-formatted to follow the PALM coding standard
28!
29! 4370 2020-01-10 14:00:44Z raasch
30! unused variables removed
31!
32! 4366 2020-01-09 08:12:43Z raasch
33! vectorized routines added
34!
35! 4182 2019-08-22 15:20:23Z scharf
36! Corrected "Former revisions" section
37!
38! 3725 2019-02-07 10:11:02Z raasch
39! GOTO statements replaced, file re-formatted corresponding to coding standards
40!
41! Revision 1.1  2003/03/12 16:41:59  raasch
42! Initial revision
43!
44!
45! Description:
46! ------------
47!> Fast Fourier transformation developed by Clive Temperton, ECMWF.
48!--------------------------------------------------------------------------------------------------!
49 MODULE temperton_fft
50
51    USE kinds
52
53    IMPLICIT NONE
54
55    PRIVATE
56
57!
58!-- No interfaces for the serial routines, because these are still writte in FORTRAN77
59    INTERFACE fft991cy_vec
60       MODULE PROCEDURE fft991cy_vec
61    END INTERFACE fft991cy_vec
62
63    INTERFACE qpassm_vec
64       MODULE PROCEDURE qpassm_vec
65    END INTERFACE qpassm_vec
66
67    INTERFACE rpassm_vec
68       MODULE PROCEDURE rpassm_vec
69    END INTERFACE rpassm_vec
70
71    PUBLIC set99, fft991cy, fft991cy_vec
72
73    INTEGER(iwp), PARAMETER ::  nfft =  32  !< maximum length of calls to *fft
74    INTEGER(iwp), PARAMETER ::  nout =   6  !< standard output stream
75
76CONTAINS
77
78!--------------------------------------------------------------------------------------------------!
79! Description:
80! ------------
81!> Calls Fortran-versions of fft's.
82!>
83!> Method:
84!>
85!> Subroutine 'fft991cy' - multiple fast real periodic transform supersedes previous routine
86!> 'fft991cy'.
87!>
88!> Real transform of length n performed by removing redundant operations from complex transform of
89!> length n.
90!>
91!> a       is the array containing input & output data.
92!> work    is an area of size (n+1)*min(lot,nfft).
93!> trigs   is a previously prepared list of trig function values.
94!> ifax    is a previously prepared list of factors of n.
95!> inc     is the increment within each data 'vector'
96!>         (e.g. inc=1 for consecutively stored data).
97!> jump    is the increment between the start of each data vector.
98!> n       is the length of the data vectors.
99!> lot     is the number of data vectors.
100!> isign = +1 for transform from spectral to gridpoint
101!>       = -1 for transform from gridpoint to spectral.
102!>
103!> ordering of coefficients:
104!> a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
105!> where b(0)=b(n/2)=0; (n+2) locations required.
106!>
107!> ordering of data:
108!> x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
109!>
110!> Vectorization is achieved on cray by doing the transforms in parallel.
111!>
112!> n must be composed of factors 2,3 & 5 but does not have to be even.
113!>
114!> definition of transforms:
115!>
116!> isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n))
117!> where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
118!>
119!> isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n))
120!> b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n))
121!>
122!> calls fortran-versions of fft's  !!!
123!> dimension a(n),work(n),trigs(n),ifax(1)
124!--------------------------------------------------------------------------------------------------!
125 SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
126
127    USE kinds
128
129    IMPLICIT NONE
130
131    INTEGER(iwp) ::  inc    !<
132    INTEGER(iwp) ::  isign  !<
133    INTEGER(iwp) ::  jump   !<
134    INTEGER(iwp) ::  lot    !<
135    INTEGER(iwp) ::  n      !<
136
137    INTEGER(iwp) ::  i        !<
138    INTEGER(iwp) ::  ia       !<
139    INTEGER(iwp) ::  ibase    !<
140    INTEGER(iwp) ::  ierr     !<
141    INTEGER(iwp) ::  ifac     !<
142    INTEGER(iwp) ::  igo      !<
143    INTEGER(iwp) ::  ii       !<
144    INTEGER(iwp) ::  ifax(*)  !<
145    INTEGER(iwp) ::  istart   !<
146    INTEGER(iwp) ::  ix       !<
147    INTEGER(iwp) ::  iz       !<
148    INTEGER(iwp) ::  j        !<
149    INTEGER(iwp) ::  jbase    !<
150    INTEGER(iwp) ::  jj       !<
151    INTEGER(iwp) ::  k        !<
152    INTEGER(iwp) ::  la       !<
153    INTEGER(iwp) ::  nb       !<
154    INTEGER(iwp) ::  nblox    !<
155    INTEGER(iwp) ::  nfax     !<
156    INTEGER(iwp) ::  nvex     !<
157    INTEGER(iwp) ::  nx       !<
158
159    REAL(wp)     ::  a(*)      !<
160    REAL(wp)     ::  trigs(*)  !<
161    REAL(wp)     ::  work(*)   !<
162
163
164    IF ( ifax(10) /= n )  CALL set99( trigs, ifax, n )
165    nfax = ifax(1)
166    nx   = n + 1
167    IF ( MOD( n, 2) == 1 )  nx = n
168    nblox = 1 + (lot-1) / nfft
169    nvex = lot - ( nblox - 1 ) * nfft
170
171
172    IF ( isign == 1 )  THEN
173!
174!--    Backward fft: spectral to gridpoint transform
175
176       istart = 1
177       DO  nb = 1, nblox
178
179          ia = istart
180          i  = istart
181          !DIR$ IVDEP
182          !CDIR NODEP
183          !OCL NOVREC
184          DO  j = 1, nvex
185             a(i+inc) = 0.5_wp * a(i)
186             i        = i + jump
187          ENDDO
188          IF ( MOD(n,2) /= 1 )  THEN
189             i = istart + n * inc
190             DO  j   = 1, nvex
191                a(i) = 0.5_wp * a(i)
192                i    = i + jump
193             ENDDO
194          ENDIF
195          ia = istart + inc
196          la = 1
197          igo = + 1
198
199          DO  k = 1, nfax
200
201             ifac = ifax(k+1)
202             ierr = -1
203             IF ( igo /= -1 )  THEN
204                CALL rpassm( a(ia), a(ia+la*inc), work(1), work(ifac*la+1), trigs, inc, 1, jump,   &
205                             nx, nvex, n, ifac, la, ierr )
206             ELSE
207                CALL rpassm( work(1), work(la+1), a(ia), a(ia+ifac*la*inc), trigs, 1, inc, nx,     &
208                             jump, nvex, n, ifac, la, ierr )
209             ENDIF
210!
211!--          Following messages shouldn't appear in PALM applications
212             IF ( ierr /= 0 )  THEN
213                SELECT CASE (ierr)
214                   CASE (:-1)
215                      WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
216                   CASE (0)
217                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
218                   CASE (1:)
219                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
220                END SELECT
221                RETURN
222             ENDIF
223
224             la  = ifac * la
225             igo = -igo
226             ia  = istart
227
228          ENDDO
229!
230!--       If necessary, copy results back to a
231          IF ( MOD( nfax, 2 ) /= 0 )  THEN
232             ibase = 1
233             jbase = ia
234             DO  jj = 1, nvex
235                i = ibase
236                j = jbase
237                DO  ii = 1, n
238                   a(j) = work(i)
239                   i    = i + 1
240                   j    = j + inc
241                ENDDO
242                ibase = ibase + nx
243                jbase = jbase + jump
244             ENDDO
245          ENDIF
246!
247!--       Fill in zeros at end
248          ix = istart + n * inc
249          !DIR$ IVDEP
250          !CDIR NODEP
251          !OCL NOVREC
252          DO  j = 1, nvex
253             a(ix)     = 0.0_wp
254             a(ix+inc) = 0.0_wp
255             ix        = ix + jump
256          ENDDO
257          istart = istart + nvex * jump
258          nvex = nfft
259
260       ENDDO
261
262    ELSEIF ( isign == -1 )  THEN
263!
264!--    Forward fft: gridpoint to spectral transform
265       istart = 1
266       DO  nb = 1, nblox
267          ia  = istart
268          la  = n
269          igo = + 1
270
271          DO  k = 1, nfax
272
273             ifac = ifax(nfax+2-k)
274             la   = la / ifac
275             ierr = -1
276             IF ( igo /= -1 )  THEN
277                CALL qpassm( a(ia), a(ia+ifac*la*inc), work(1), work(la+1), trigs, inc, 1, jump,   &
278                             nx, nvex, n, ifac, la, ierr )
279             ELSE
280                CALL qpassm( work(1), work(ifac*la+1), a(ia), a(ia+la*inc), trigs, 1, inc, nx,     &
281                             jump, nvex, n, ifac, la, ierr )
282             ENDIF
283!
284!--          Following messages shouldn't appear in PALM applications
285             IF ( ierr /= 0 )  THEN
286                SELECT CASE (ierr)
287                   CASE (:-1)
288                      WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
289                   CASE (0)
290                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
291                   CASE (1:)
292                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
293                END SELECT
294                RETURN
295             ENDIF
296
297             igo = -igo
298             ia = istart + inc
299
300          ENDDO
301!
302!--       If necessary, copy results back to a
303          IF ( MOD(nfax,2) /= 0 )  THEN
304             ibase = 1
305             jbase = ia
306             DO  jj = 1, nvex
307                i = ibase
308                j = jbase
309                DO  ii = 1, n
310                   a(j) = work(i)
311                   i    = i + 1
312                   j    = j + inc
313                ENDDO
314                ibase = ibase + nx
315                jbase = jbase + jump
316             ENDDO
317          ENDIF
318!
319!--       Shift a(0) and fill in zero imag parts
320          ix = istart
321          !DIR$ IVDEP
322          !CDIR NODEP
323          !OCL NOVREC
324          DO  j = 1, nvex
325             a(ix)     = a(ix+inc)
326             a(ix+inc) = 0.0_wp
327             ix        = ix + jump
328          ENDDO
329          IF ( MOD(n,2) /= 1 )  THEN
330             iz = istart + (n+1) * inc
331             DO  j = 1, nvex
332                a(iz) = 0.0_wp
333                iz    = iz + jump
334             ENDDO
335          ENDIF
336          istart = istart + nvex * jump
337          nvex = nfft
338
339       ENDDO
340
341    ENDIF
342
343 END SUBROUTINE fft991cy
344
345!--------------------------------------------------------------------------------------------------!
346! Description:
347! ------------
348!> Performs one pass through data as part of multiple real fft (fourier analysis) routine.
349!>
350!> Method:
351!>
352!> a       is first real input vector
353!>         equivalence b(1) with a(ifac*la*inc1+1)
354!> c       is first real output vector
355!>         equivalence d(1) with c(la*inc2+1)
356!> trigs   is a precalculated list of sines & cosines
357!> inc1    is the addressing increment for a
358!> inc2    is the addressing increment for c
359!> inc3    is the increment between input vectors a
360!> inc4    is the increment between output vectors c
361!> lot     is the number of vectors
362!> n       is the length of the vectors
363!> ifac    is the current factor of n
364!>         la = n/(product of factors used so far)
365!> ierr    is an error indicator:
366!>         0 - pass completed without error
367!>         1 - lot greater than nfft
368!>         2 - ifac not catered for
369!>         3 - ifac only catered for if la=n/ifac
370!--------------------------------------------------------------------------------------------------!
371 SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
372
373    USE kinds
374
375    IMPLICIT NONE
376
377    INTEGER(iwp) ::  ierr  !<
378    INTEGER(iwp) ::  ifac  !<
379    INTEGER(iwp) ::  inc1  !<
380    INTEGER(iwp) ::  inc2  !<
381    INTEGER(iwp) ::  inc3  !<
382    INTEGER(iwp) ::  inc4  !<
383    INTEGER(iwp) ::  la    !<
384    INTEGER(iwp) ::  lot   !<
385    INTEGER(iwp) ::  n     !<
386
387    INTEGER(iwp) ::  i      !<
388    INTEGER(iwp) ::  ia     !<
389    INTEGER(iwp) ::  ib     !<
390    INTEGER(iwp) ::  ibase  !<
391    INTEGER(iwp) ::  ic     !<
392    INTEGER(iwp) ::  id     !<
393    INTEGER(iwp) ::  ie     !<
394    INTEGER(iwp) ::  if     !<
395    INTEGER(iwp) ::  ig     !<
396    INTEGER(iwp) ::  igo    !<
397    INTEGER(iwp) ::  ih     !<
398    INTEGER(iwp) ::  iink   !<
399    INTEGER(iwp) ::  ijk    !<
400    INTEGER(iwp) ::  ijump  !<
401    INTEGER(iwp) ::  j      !<
402    INTEGER(iwp) ::  ja     !<
403    INTEGER(iwp) ::  jb     !<
404    INTEGER(iwp) ::  jbase  !<
405    INTEGER(iwp) ::  jc     !<
406    INTEGER(iwp) ::  jd     !<
407    INTEGER(iwp) ::  je     !<
408    INTEGER(iwp) ::  jf     !<
409    INTEGER(iwp) ::  jink   !<
410    INTEGER(iwp) ::  k      !<
411    INTEGER(iwp) ::  kb     !<
412    INTEGER(iwp) ::  kc     !<
413    INTEGER(iwp) ::  kd     !<
414    INTEGER(iwp) ::  ke     !<
415    INTEGER(iwp) ::  kf     !<
416    INTEGER(iwp) ::  kstop  !<
417    INTEGER(iwp) ::  l      !<
418    INTEGER(iwp) ::  m      !<
419
420!
421!-- Arrays are dimensioned with n
422    REAL(wp) ::  a(*)      !<
423    REAL(wp) ::  b(*)      !<
424    REAL(wp) ::  c(*)      !<
425    REAL(wp) ::  d(*)      !<
426    REAL(wp) ::  trigs(*)  !<
427
428    REAL(wp) ::  a0      !<
429    REAL(wp) ::  a1      !<
430    REAL(wp) ::  a10     !<
431    REAL(wp) ::  a11     !<
432    REAL(wp) ::  a2      !<
433    REAL(wp) ::  a20     !<
434    REAL(wp) ::  a21     !<
435    REAL(wp) ::  a3      !<
436    REAL(wp) ::  a4      !<
437    REAL(wp) ::  a5      !<
438    REAL(wp) ::  a6      !<
439    REAL(wp) ::  b0      !<
440    REAL(wp) ::  b1      !<
441    REAL(wp) ::  b10     !<
442    REAL(wp) ::  b11     !<
443    REAL(wp) ::  b2      !<
444    REAL(wp) ::  b20     !<
445    REAL(wp) ::  b21     !<
446    REAL(wp) ::  b3      !<
447    REAL(wp) ::  b4      !<
448    REAL(wp) ::  b5      !<
449    REAL(wp) ::  b6      !<
450    REAL(wp) ::  c1      !<
451    REAL(wp) ::  c2      !<
452    REAL(wp) ::  c3      !<
453    REAL(wp) ::  c4      !<
454    REAL(wp) ::  c5      !<
455    REAL(wp) ::  qrt5    !<
456    REAL(wp) ::  s1      !<
457    REAL(wp) ::  s2      !<
458    REAL(wp) ::  s3      !<
459    REAL(wp) ::  s4      !<
460    REAL(wp) ::  s5      !<
461    REAL(wp) ::  sin36   !<
462    REAL(wp) ::  sin45   !<
463    REAL(wp) ::  sin60   !<
464    REAL(wp) ::  sin72   !<
465    REAL(wp) ::  z       !<
466    REAL(wp) ::  zqrt5   !<
467    REAL(wp) ::  zsin36  !<
468    REAL(wp) ::  zsin45  !<
469    REAL(wp) ::  zsin60  !<
470    REAL(wp) ::  zsin72  !<
471
472
473    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
474          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
475
476
477    ierr = 0
478
479    IF ( lot > nfft )  THEN
480       ierr = 1
481       RETURN
482    ENDIF
483
484    m = n / ifac
485    iink  = la * inc1
486    jink  = la * inc2
487    ijump = (ifac - 1) * iink
488    kstop = ( n - ifac ) / ( 2 * ifac )
489
490    ibase = 0
491    jbase = 0
492    igo = ifac - 1
493    IF ( igo == 7 )  igo = 6
494
495    IF (igo < 1 .OR. igo > 6 )  THEN
496       ierr = 2
497       RETURN
498    ENDIF
499
500
501    SELECT CASE ( igo )
502
503!
504!--    Coding for factor 2
505       CASE ( 1 )
506          ia = 1
507          ib = ia + iink
508          ja = 1
509          jb = ja + ( 2 * m - la ) * inc2
510
511          IF ( la /= m )  THEN
512
513             DO  l = 1, la
514                i = ibase
515                j = jbase
516                !DIR$ IVDEP
517                !CDIR NODEP
518                !OCL NOVREC
519                DO  ijk = 1, lot
520                   c(ja+j) = a(ia+i) + a(ib+i)
521                   c(jb+j) = a(ia+i) - a(ib+i)
522                   i = i + inc3
523                   j = j + inc4
524                ENDDO
525                ibase = ibase + inc1
526                jbase = jbase + inc2
527             ENDDO
528
529             ja    = ja + jink
530             jink  = 2 * jink
531             jb    = jb - jink
532             ibase = ibase + ijump
533             ijump = 2 * ijump + iink
534
535             IF ( ja /= jb )  THEN
536
537                DO  k = la, kstop, la
538                   kb = k + k
539                   c1 = trigs(kb+1)
540                   s1 = trigs(kb+2)
541                   jbase = 0
542                   DO  l = 1, la
543                      i = ibase
544                      j = jbase
545                      !DIR$ IVDEP
546                      !CDIR NODEP
547                      !OCL NOVREC
548                      DO  ijk = 1, lot
549                         c(ja+j) = a(ia+i) + ( c1 * a(ib+i) + s1 * b(ib+i) )
550                         c(jb+j) = a(ia+i) - ( c1 * a(ib+i) + s1 * b(ib+i) )
551                         d(ja+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) + b(ia+i)
552                         d(jb+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) - b(ia+i)
553                         i = i + inc3
554                         j = j + inc4
555                      ENDDO
556                      ibase = ibase + inc1
557                      jbase = jbase + inc2
558                   ENDDO
559
560                   ibase = ibase + ijump
561                   ja = ja + jink
562                   jb = jb - jink
563                ENDDO
564
565                IF ( ja > jb )  RETURN
566
567             ENDIF
568
569             jbase = 0
570             DO  l = 1, la
571                i = ibase
572                j = jbase
573                !DIR$ IVDEP
574                !CDIR NODEP
575                !OCL NOVREC
576                DO  ijk = 1, lot
577                   c(ja+j) = a(ia+i)
578                   d(ja+j) = -a(ib+i)
579                   i = i + inc3
580                   j = j + inc4
581                ENDDO
582                ibase = ibase + inc1
583                jbase = jbase + inc2
584             ENDDO
585
586          ELSE
587
588             z = 1.0_wp / REAL( n, KIND = wp )
589             DO  l = 1, la
590                i = ibase
591                j = jbase
592                !DIR$ IVDEP
593                !CDIR NODEP
594                !OCL NOVREC
595                DO  ijk = 1, lot
596                   c(ja+j) = z * ( a(ia+i) + a(ib+i) )
597                   c(jb+j) = z * ( a(ia+i) - a(ib+i) )
598                   i = i + inc3
599                   j = j + inc4
600                ENDDO
601                ibase = ibase + inc1
602                jbase = jbase + inc2
603             ENDDO
604
605          ENDIF
606
607!
608!--    Coding for factor 3
609       CASE ( 2 )
610
611          ia = 1
612          ib = ia + iink
613          ic = ib + iink
614          ja = 1
615          jb = ja + ( 2 * m - la) * inc2
616          jc = jb
617
618          IF ( la /= m )  THEN
619
620             DO  l = 1, la
621                i = ibase
622                j = jbase
623                !DIR$ IVDEP
624                !CDIR NODEP
625                !OCL NOVREC
626                DO  ijk = 1, lot
627                   c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
628                   c(jb+j) = a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) )
629                   d(jb+j) = sin60 * ( a(ic+i) - a(ib+i) )
630                   i = i + inc3
631                   j = j + inc4
632                ENDDO
633                ibase = ibase + inc1
634                jbase = jbase + inc2
635             ENDDO
636
637             ja = ja + jink
638             jink = 2 * jink
639             jb = jb + jink
640             jc = jc - jink
641             ibase = ibase + ijump
642             ijump = 2 * ijump + iink
643
644             IF ( ja /= jc )  THEN
645
646                DO  k = la, kstop, la
647                   kb = k + k
648                   kc = kb + kb
649                   c1 = trigs(kb+1)
650                   s1 = trigs(kb+2)
651                   c2 = trigs(kc+1)
652                   s2 = trigs(kc+2)
653
654                   jbase = 0
655                   DO  l = 1, la
656                      i = ibase
657                      j = jbase
658                      !DIR$ IVDEP
659                      !CDIR NODEP
660                      !OCL NOVREC
661                      DO  ijk = 1, lot
662                         a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c2 * a(ic+i) + s2 * b(ic+i) )
663                         b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c2 * b(ic+i) - s2 * a(ic+i) )
664                         a2 = a(ia+i) - 0.5_wp * a1
665                         b2 = b(ia+i) - 0.5_wp * b1
666                         a3 = sin60 * ( ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c2 * a(ic+i) + s2      &
667                              * b(ic+i) ) )
668                         b3 = sin60 * ( ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c2 * b(ic+i) - s2      &
669                              * a(ic+i) ) )
670                         c(ja+j) = a(ia+i) + a1
671                         d(ja+j) = b(ia+i) + b1
672                         c(jb+j) = a2 + b3
673                         d(jb+j) = b2 - a3
674                         c(jc+j) = a2 - b3
675                         d(jc+j) = - ( b2 + a3 )
676                         i = i + inc3
677                         j = j + inc4
678                      ENDDO
679                      ibase = ibase + inc1
680                      jbase = jbase + inc2
681                   ENDDO
682                   ibase = ibase + ijump
683                   ja = ja + jink
684                   jb = jb + jink
685                   jc = jc - jink
686                ENDDO
687
688                IF ( ja > jc )  RETURN
689
690             ENDIF
691
692             jbase = 0
693             DO  l = 1, la
694                i = ibase
695                j = jbase
696                !DIR$ IVDEP
697                !CDIR NODEP
698                !OCL NOVREC
699                DO  ijk = 1, lot
700                   c(ja+j) = a(ia+i) + 0.5_wp * ( a(ib+i) - a(ic+i) )
701                   d(ja+j) = -sin60 * ( a(ib+i) + a(ic+i) )
702                   c(jb+j) = a(ia+i) - ( a(ib+i) - a(ic+i) )
703                   i = i + inc3
704                   j = j + inc4
705                ENDDO
706                ibase = ibase + inc1
707                jbase = jbase + inc2
708             ENDDO
709
710          ELSE
711
712             z = 1.0_wp / REAL( n, KIND = wp )
713             zsin60 = z * sin60
714             DO  l = 1, la
715                i = ibase
716                j = jbase
717                !DIR$ IVDEP
718                !CDIR NODEP
719                !OCL NOVREC
720                DO  ijk = 1, lot
721                   c(ja+j) = z * ( a(ia+i) + ( a (ib+i) + a(ic+i) ) )
722                   c(jb+j) = z * ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )
723                   d(jb+j) = zsin60 * ( a(ic+i) -a(ib+i) )
724                   i = i + inc3
725                   j = j + inc4
726                ENDDO
727                ibase = ibase + inc1
728                jbase = jbase + inc2
729             ENDDO
730
731          ENDIF
732
733!
734!--    Coding for factor 4
735       CASE ( 3 )
736
737          ia = 1
738          ib = ia + iink
739          ic = ib + iink
740          id = ic + iink
741          ja = 1
742          jb = ja + ( 2 * m - la ) * inc2
743          jc = jb + 2 * m * inc2
744          jd = jb
745
746          IF ( la /= m )  THEN
747
748             DO  l = 1, la
749
750                i = ibase
751                j = jbase
752                !DIR$ IVDEP
753                !CDIR NODEP
754                !OCL NOVREC
755                DO  ijk = 1, lot
756                   c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) )
757                   c(jc+j) = ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) )
758                   c(jb+j) = a(ia+i) - a(ic+i)
759                   d(jb+j) = a(id+i) - a(ib+i)
760                   i = i + inc3
761                   j = j + inc4
762                ENDDO
763                ibase = ibase + inc1
764                jbase = jbase + inc2
765             ENDDO
766
767             ja    = ja + jink
768             jink  = 2 * jink
769             jb    = jb + jink
770             jc    = jc - jink
771             jd    = jd - jink
772             ibase = ibase + ijump
773             ijump = 2 * ijump + iink
774
775             IF ( jb /= jc )  THEN
776
777                DO  k = la, kstop, la
778                   kb = k + k
779                   kc = kb + kb
780                   kd = kc + kb
781                   c1 = trigs(kb+1)
782                   s1 = trigs(kb+2)
783                   c2 = trigs(kc+1)
784                   s2 = trigs(kc+2)
785                   c3 = trigs(kd+1)
786                   s3 = trigs(kd+2)
787                   jbase = 0
788                   DO  l = 1, la
789                      i = ibase
790                      j = jbase
791                      !DIR$ IVDEP
792                      !CDIR NODEP
793                      !OCL NOVREC
794                      DO  ijk = 1, lot
795                         a0 = a(ia+i) + ( c2 * a(ic+i) + s2 * b(ic+i) )
796                         a2 = a(ia+i) - ( c2 * a(ic+i) + s2 * b(ic+i) )
797                         a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c3 * a(id+i) + s3 * b(id+i) )
798                         a3 = ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c3 * a(id+i) + s3 * b(id+i) )
799                         b0 = b(ia+i) + ( c2 * b(ic+i) - s2 * a(ic+i) )
800                         b2 = b(ia+i) - ( c2 * b(ic+i) - s2 * a(ic+i) )
801                         b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c3 * b(id+i) - s3 * a(id+i) )
802                         b3 = ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c3 * b(id+i) - s3 * a(id+i) )
803                         c(ja+j) = a0 + a1
804                         c(jc+j) = a0 - a1
805                         d(ja+j) = b0 + b1
806                         d(jc+j) = b1 - b0
807                         c(jb+j) = a2 + b3
808                         c(jd+j) = a2 - b3
809                         d(jb+j) = b2 - a3
810                         d(jd+j) = - ( b2 + a3 )
811                         i = i + inc3
812                         j = j + inc4
813                      ENDDO
814                      ibase = ibase + inc1
815                      jbase = jbase + inc2
816                   ENDDO
817                   ibase = ibase + ijump
818                   ja = ja + jink
819                   jb = jb + jink
820                   jc = jc - jink
821                   jd = jd - jink
822                ENDDO
823
824                IF ( jb > jc )  RETURN
825
826             ENDIF
827
828             sin45 = SQRT( 0.5_wp )
829             jbase = 0
830             DO  l = 1, la
831                i = ibase
832                j = jbase
833                !DIR$ IVDEP
834                !CDIR NODEP
835                !OCL NOVREC
836                DO  ijk = 1, lot
837                   c(ja+j) = a(ia+i) + sin45 * ( a(ib+i) - a(id+i) )
838                   c(jb+j) = a(ia+i) - sin45 * ( a(ib+i) - a(id+i) )
839                   d(ja+j) = - a(ic+i) - sin45 * ( a(ib+i) + a(id+i) )
840                   d(jb+j) = a(ic+i) - sin45 * ( a(ib+i) + a(id+i) )
841                   i = i + inc3
842                   j = j + inc4
843                ENDDO
844                ibase = ibase + inc1
845                jbase = jbase + inc2
846             ENDDO
847
848          ELSE
849
850             z = 1.0_wp / REAL( n, KIND = wp )
851             DO  l = 1, la
852                i = ibase
853                j = jbase
854                !DIR$ IVDEP
855                !CDIR NODEP
856                !OCL NOVREC
857                DO  ijk = 1, lot
858                   c(ja+j) = z * ( ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) )
859                   c(jc+j) = z * ( ( a(ia+i) + a(ic+i) ) - (a (ib+i) + a(id+i) ) )
860                   c(jb+j) = z * ( a(ia+i) - a(ic+i) )
861                   d(jb+j) = z * ( a(id+i) - a(ib+i) )
862                   i = i + inc3
863                   j = j + inc4
864                ENDDO
865                ibase = ibase + inc1
866                jbase = jbase + inc2
867             ENDDO
868
869          ENDIF
870
871!
872!--    Coding for factor 5
873       CASE ( 4 )
874
875          ia = 1
876          ib = ia + iink
877          ic = ib + iink
878          id = ic + iink
879          ie = id + iink
880          ja = 1
881          jb = ja + ( 2 * m - la ) * inc2
882          jc = jb + 2 * m * inc2
883          jd = jc
884          je = jb
885
886          IF ( la /= m )  THEN
887
888             DO  l = 1, la
889                i = ibase
890                j = jbase
891                !DIR$ IVDEP
892                !CDIR NODEP
893                !OCL NOVREC
894                DO  ijk = 1, lot
895                   a1 = a(ib+i) + a(ie+i)
896                   a3 = a(ib+i) - a(ie+i)
897                   a2 = a(ic+i) + a(id+i)
898                   a4 = a(ic+i) - a(id+i)
899                   a5 = a(ia+i) - 0.25_wp * ( a1 + a2 )
900                   a6 = qrt5 * ( a1 - a2 )
901                   c(ja+j) = a(ia+i) + ( a1 + a2 )
902                   c(jb+j) = a5 + a6
903                   c(jc+j) = a5 - a6
904                   d(jb+j) = - sin72 * a3 - sin36 * a4
905                   d(jc+j) = - sin36 * a3 + sin72 * a4
906                   i = i + inc3
907                   j = j + inc4
908                ENDDO
909                ibase = ibase + inc1
910                jbase = jbase + inc2
911             ENDDO
912
913             ja = ja + jink
914             jink = 2 * jink
915             jb = jb + jink
916             jc = jc + jink
917             jd = jd - jink
918             je = je - jink
919             ibase = ibase + ijump
920             ijump = 2 * ijump + iink
921
922             IF ( jb /= jd )  THEN
923
924                DO  k = la, kstop, la
925                   kb = k + k
926                   kc = kb + kb
927                   kd = kc + kb
928                   ke = kd + kb
929                   c1 = trigs(kb+1)
930                   s1 = trigs(kb+2)
931                   c2 = trigs(kc+1)
932                   s2 = trigs(kc+2)
933                   c3 = trigs(kd+1)
934                   s3 = trigs(kd+2)
935                   c4 = trigs(ke+1)
936                   s4 = trigs(ke+2)
937                   jbase = 0
938                   DO  l = 1, la
939                      i = ibase
940                      j = jbase
941                      !DIR$ IVDEP
942                      !CDIR NODEP
943                      !OCL NOVREC
944                      DO  ijk = 1, lot
945                         a1  = (c1 * a(ib+i) + s1 * b(ib+i) ) + ( c4 * a(ie+i) + s4 * b(ie+i) )
946                         a3  = (c1 * a(ib+i) + s1 * b(ib+i) ) - ( c4 * a(ie+i) + s4 * b(ie+i) )
947                         a2  = (c2 * a(ic+i) + s2 * b(ic+i) ) + ( c3 * a(id+i) + s3 * b(id+i) )
948                         a4  = (c2 * a(ic+i) + s2 * b(ic+i) ) - ( c3 * a(id+i) + s3 * b(id+i) )
949                         b1  = (c1 * b(ib+i) - s1 * a(ib+i) ) + ( c4 * b(ie+i) - s4 * a(ie+i) )
950                         b3  = (c1 * b(ib+i) - s1 * a(ib+i) ) - ( c4 * b(ie+i) - s4 * a(ie+i) )
951                         b2  = (c2 * b(ic+i) - s2 * a(ic+i) ) + ( c3 * b(id+i) - s3 * a(id+i) )
952                         b4  = (c2 * b(ic+i) - s2 * a(ic+i) ) - ( c3 * b(id+i) - s3 * a(id+i) )
953                         a5  = a(ia+i) - 0.25_wp * ( a1 + a2 )
954                         a6  = qrt5 * ( a1 - a2 )
955                         b5  = b(ia+i) - 0.25_wp * ( b1 + b2 )
956                         b6  = qrt5 * ( b1 - b2 )
957                         a10 = a5 + a6
958                         a20 = a5 - a6
959                         b10 = b5 + b6
960                         b20 = b5 - b6
961                         a11 = sin72 * b3 + sin36 * b4
962                         a21 = sin36 * b3 - sin72 * b4
963                         b11 = sin72 * a3 + sin36 * a4
964                         b21 = sin36 * a3 - sin72 * a4
965                         c(ja+j) = a(ia+i) + ( a1 + a2 )
966                         c(jb+j) = a10 + a11
967                         c(je+j) = a10 - a11
968                         c(jc+j) = a20 + a21
969                         c(jd+j) = a20 - a21
970                         d(ja+j) = b(ia+i) + ( b1 + b2 )
971                         d(jb+j) = b10 - b11
972                         d(je+j) = - ( b10 + b11 )
973                         d(jc+j) = b20 - b21
974                         d(jd+j) = - ( b20 + b21 )
975                         i = i + inc3
976                         j = j + inc4
977                      ENDDO
978                      ibase = ibase + inc1
979                      jbase = jbase + inc2
980                   ENDDO
981                   ibase = ibase + ijump
982                   ja = ja + jink
983                   jb = jb + jink
984                   jc = jc + jink
985                   jd = jd - jink
986                   je = je - jink
987                ENDDO
988
989                IF ( jb > jd )  RETURN
990
991             ENDIF
992
993             jbase = 0
994             DO  l = 1, la
995                i = ibase
996                j = jbase
997                !DIR$ IVDEP
998                !CDIR NODEP
999                !OCL NOVREC
1000                DO  ijk = 1, lot
1001                   a1 = a(ib+i) + a(ie+i)
1002                   a3 = a(ib+i) - a(ie+i)
1003                   a2 = a(ic+i) + a(id+i)
1004                   a4 = a(ic+i) - a(id+i)
1005                   a5 = a(ia+i) + 0.25_wp * ( a3 - a4 )
1006                   a6 = qrt5 * ( a3 + a4 )
1007                   c(ja+j) = a5 + a6
1008                   c(jb+j) = a5 - a6
1009                   c(jc+j) = a(ia+i) - ( a3 - a4 )
1010                   d(ja+j) = - sin36 * a1 - sin72 * a2
1011                   d(jb+j) = - sin72 * a1 + sin36 * a2
1012                   i = i + inc3
1013                   j = j + inc4
1014                ENDDO
1015                ibase = ibase + inc1
1016                jbase = jbase + inc2
1017             ENDDO
1018
1019          ELSE
1020
1021             z = 1.0_wp / REAL( n, KIND = wp )
1022             zqrt5  = z * qrt5
1023             zsin36 = z * sin36
1024             zsin72 = z * sin72
1025             DO  l = 1, la
1026                i = ibase
1027                j = jbase
1028                !DIR$ IVDEP
1029                !CDIR NODEP
1030                !OCL NOVREC
1031                DO  ijk = 1, lot
1032                   a1 = a(ib+i) + a(ie+i)
1033                   a3 = a(ib+i) - a(ie+i)
1034                   a2 = a(ic+i) + a(id+i)
1035                   a4 = a(ic+i) - a(id+i)
1036                   a5 = z * ( a (ia+i) - 0.25_wp * ( a1 + a2 ) )
1037                   a6 = zqrt5 * ( a1 - a2 )
1038                   c(ja+j) = z * ( a(ia+i) + (a1+a2) )
1039                   c(jb+j) = a5 + a6
1040                   c(jc+j) = a5 - a6
1041                   d(jb+j) = - zsin72 * a3 - zsin36 * a4
1042                   d(jc+j) = - zsin36 * a3 + zsin72 * a4
1043                   i = i + inc3
1044                   j = j + inc4
1045                ENDDO
1046                ibase = ibase + inc1
1047                jbase = jbase + inc2
1048             ENDDO
1049
1050          ENDIF
1051
1052!
1053!--    Coding for factor 6
1054       CASE ( 5 )
1055
1056          ia = 1
1057          ib = ia + iink
1058          ic = ib + iink
1059          id = ic + iink
1060          ie = id + iink
1061          if = ie + iink
1062          ja = 1
1063          jb = ja + ( 2 * m - la ) * inc2
1064          jc = jb + 2 * m * inc2
1065          jd = jc + 2 * m * inc2
1066          je = jc
1067          jf = jb
1068
1069          IF ( la /= m )  THEN
1070
1071             DO  l = 1, la
1072                i = ibase
1073                j = jbase
1074                !DIR$ IVDEP
1075                !CDIR NODEP
1076                !OCL NOVREC
1077                DO  ijk = 1, lot
1078                   a11 = ( a (ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) )
1079                   c(ja+j) = ( a(ia+i) + a(id+i) ) + a11
1080                   c(jc+j) = ( a(ia+i) + a(id+i) - 0.5_wp * a11 )
1081                   d(jc+j) = sin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) )
1082                   a11 = ( a (ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) )
1083                   c(jb+j) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11
1084                   d(jb+j) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
1085                   c(jd+j) = ( a (ia+i) - a(id+i) ) + a11
1086                   i = i + inc3
1087                   j = j + inc4
1088                ENDDO
1089                ibase = ibase + inc1
1090                jbase = jbase + inc2
1091             ENDDO
1092             ja    = ja + jink
1093             jink  = 2 * jink
1094             jb    = jb + jink
1095             jc    = jc + jink
1096             jd    = jd - jink
1097             je    = je - jink
1098             jf    = jf - jink
1099             ibase = ibase + ijump
1100             ijump = 2 * ijump + iink
1101
1102             IF ( jc /= jd )  THEN
1103
1104                DO  k = la, kstop, la
1105                   kb = k + k
1106                   kc = kb + kb
1107                   kd = kc + kb
1108                   ke = kd + kb
1109                   kf = ke + kb
1110                   c1 = trigs(kb+1)
1111                   s1 = trigs(kb+2)
1112                   c2 = trigs(kc+1)
1113                   s2 = trigs(kc+2)
1114                   c3 = trigs(kd+1)
1115                   s3 = trigs(kd+2)
1116                   c4 = trigs(ke+1)
1117                   s4 = trigs(ke+2)
1118                   c5 = trigs(kf+1)
1119                   s5 = trigs(kf+2)
1120                   jbase = 0
1121                   DO  l = 1, la
1122                      i = ibase
1123                      j = jbase
1124                      !DIR$ IVDEP
1125                      !CDIR NODEP
1126                      !OCL NOVREC
1127                      DO  ijk = 1, lot
1128                         a1  = c1 * a(ib+i) + s1 * b(ib+i)
1129                         b1  = c1 * b(ib+i) - s1 * a(ib+i)
1130                         a2  = c2 * a(ic+i) + s2 * b(ic+i)
1131                         b2  = c2 * b(ic+i) - s2 * a(ic+i)
1132                         a3  = c3 * a(id+i) + s3 * b(id+i)
1133                         b3  = c3 * b(id+i) - s3 * a(id+i)
1134                         a4  = c4 * a(ie+i) + s4 * b(ie+i)
1135                         b4  = c4 * b(ie+i) - s4 * a(ie+i)
1136                         a5  = c5 * a(if+i) + s5 * b(if+i)
1137                         b5  = c5 * b(if+i) - s5 * a(if+i)
1138                         a11 = ( a2 + a5 ) + ( a1 + a4 )
1139                         a20 = ( a(ia+i) + a3 ) - 0.5_wp * a11
1140                         a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) )
1141                         b11 = ( b2 + b5 ) + ( b1 + b4 )
1142                         b20 = ( b(ia+i) + b3 ) - 0.5_wp * b11
1143                         b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) )
1144                         c(ja+j) = ( a(ia+i) + a3 ) + a11
1145                         d(ja+j) = ( b(ia+i) + b3 ) + b11
1146                         c(jc+j) = a20 - b21
1147                         d(jc+j) = a21 + b20
1148                         c(je+j) = a20 + b21
1149                         d(je+j) = a21 - b20
1150                         a11 = ( a2 - a5 ) + ( a4 - a1 )
1151                         a20 = ( a(ia+i) - a3 ) - 0.5_wp * a11
1152                         a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) )
1153                         b11 = ( b5 - b2 ) - ( b4 - b1 )
1154                         b20 = ( b3 - b(ia+i) ) - 0.5_wp * b11
1155                         b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) )
1156                         c(jb+j) = a20 - b21
1157                         d(jb+j) = a21 - b20
1158                         c(jd+j) = a11 + ( a(ia+i) - a3 )
1159                         d(jd+j) = b11 + ( b3 - b(ia+i) )
1160                         c(jf+j) = a20 + b21
1161                         d(jf+j) = a21 + b20
1162                         i = i + inc3
1163                         j = j + inc4
1164                      ENDDO
1165                      ibase = ibase + inc1
1166                      jbase = jbase + inc2
1167                   ENDDO
1168                   ibase = ibase + ijump
1169                   ja = ja + jink
1170                   jb = jb + jink
1171                   jc = jc + jink
1172                   jd = jd - jink
1173                   je = je - jink
1174                   jf = jf - jink
1175                ENDDO
1176
1177                IF ( jc > jd )  RETURN
1178
1179             ENDIF
1180
1181             jbase = 0
1182             DO  l = 1, la
1183                i = ibase
1184                j = jbase
1185                !DIR$ IVDEP
1186                !CDIR NODEP
1187                !OCL NOVREC
1188                DO  ijk = 1, lot
1189                   c(ja+j) = ( a (ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) + sin60 *               &
1190                             ( a(ib+i) - a(if+i) )
1191                   d(ja+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) - sin60 * ( a(ic+i) +  &
1192                             a(ie+i) )
1193                   c(jb+j) = a(ia+i) - ( a(ic+i) - a(ie+i) )
1194                   d(jb+j) = a(id+i) - ( a(ib+i) + a(if+i) )
1195                   c(jc+j) = ( a(ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) - sin60 * ( a(ib+i) -    &
1196                             a(if+i) )
1197                   d(jc+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) + sin60 * ( a(ic+i) +  &
1198                             a(ie+i) )
1199                   i = i + inc3
1200                   j = j + inc4
1201                ENDDO
1202                ibase = ibase + inc1
1203                jbase = jbase + inc2
1204             ENDDO
1205
1206          ELSE
1207
1208             z = 1.0_wp / REAL( n, KIND = wp )
1209             zsin60 = z * sin60
1210             DO  l = 1, la
1211                i = ibase
1212                j = jbase
1213                !DIR$ IVDEP
1214                !CDIR NODEP
1215                !OCL NOVREC
1216                DO  ijk = 1, lot
1217                   a11 = ( a(ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) )
1218                   c(ja+j) = z * ( ( a(ia+i) + a(id+i) ) + a11 )
1219                   c(jc+j) = z * ( ( a(ia+i) + a(id+i) ) - 0.5_wp * a11 )
1220                   d(jc+j) = zsin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) )
1221                   a11 = ( a(ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) )
1222                   c(jb+j) = z * ( ( a(ia+i) - a(id+i) ) - 0.5_wp * a11 )
1223                   d(jb+j) = zsin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
1224                   c(jd+j) = z * ( ( a(ia+i) - a(id+i) ) + a11 )
1225                   i = i + inc3
1226                   j = j + inc4
1227                ENDDO
1228                ibase = ibase + inc1
1229                jbase = jbase + inc2
1230             ENDDO
1231
1232          ENDIF
1233
1234!
1235!--    Coding for factor 8
1236       CASE ( 6 )
1237
1238          IF ( la /= m )  THEN
1239             ierr = 3
1240             RETURN
1241          ENDIF
1242
1243          ia = 1
1244          ib = ia + iink
1245          ic = ib + iink
1246          id = ic + iink
1247          ie = id + iink
1248          if = ie + iink
1249          ig = if + iink
1250          ih = ig + iink
1251          ja = 1
1252          jb = ja + la * inc2
1253          jc = jb + 2 * m * inc2
1254          jd = jc + 2 * m * inc2
1255          je = jd + 2 * m * inc2
1256          z = 1.0_wp / REAL( n, KIND = wp )
1257          zsin45 = z * SQRT( 0.5_wp )
1258
1259          DO  l = 1, la
1260             i = ibase
1261             j = jbase
1262
1263             !DIR$ IVDEP
1264             !CDIR NODEP
1265             !OCL NOVREC
1266             DO  ijk = 1, lot
1267                c(ja+j) = z * ( ( ( a(ia+i) + a(ie+i) ) + ( a(ic+i) + a(ig+i) ) ) + ( ( a(id+i)    &
1268                          + a(ih+i) ) + ( a(ib+i) + a(if+i) ) ) )
1269                c(je+j) = z * ( ( ( a(ia+i) + a(ie+i) ) + (a(ic+i) + a(ig+i) ) ) - ( ( a(id+i)     &
1270                          + a(ih+i) ) + ( a(ib+i) + a(if+i) ) ) )
1271                c(jc+j) = z * ( ( a(ia+i) + a(ie+i) ) - ( a(ic+i) + a(ig+i) ) )
1272                d(jc+j) = z * ( ( a(id+i) + a(ih+i) ) - ( a(ib+i) + a(if+i) ) )
1273                c(jb+j) = z * ( a(ia+i) - a(ie+i) ) + zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) &
1274                          - a(ib+i) ) )
1275                c(jd+j) = z * ( a(ia+i) - a(ie+i) ) - zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) &
1276                          - a(ib+i) ) )
1277                d(jb+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) + z *         &
1278                          ( a(ig+i) - a(ic+i) )
1279                d(jd+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) - z *         &
1280                          ( a(ig+i) - a(ic+i) )
1281                i = i + inc3
1282                j = j + inc4
1283             ENDDO
1284             ibase = ibase + inc1
1285             jbase = jbase + inc2
1286          ENDDO
1287
1288    END SELECT
1289
1290 END SUBROUTINE qpassm
1291
1292!--------------------------------------------------------------------------------------------------!
1293! Description:
1294! ------------
1295!> Same as qpassm, but for backward fft
1296!--------------------------------------------------------------------------------------------------!
1297 SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
1298
1299    USE kinds
1300
1301    IMPLICIT NONE
1302
1303    INTEGER(iwp) ::  ierr  !<
1304    INTEGER(iwp) ::  ifac  !<
1305    INTEGER(iwp) ::  inc1  !<
1306    INTEGER(iwp) ::  inc2  !<
1307    INTEGER(iwp) ::  inc3  !<
1308    INTEGER(iwp) ::  inc4  !<
1309    INTEGER(iwp) ::  la    !<
1310    INTEGER(iwp) ::  lot   !<
1311    INTEGER(iwp) ::  n     !<
1312!
1313!-- Arrays are dimensioned with n
1314    INTEGER(iwp) ::  i      !<
1315    INTEGER(iwp) ::  ia     !<
1316    INTEGER(iwp) ::  ib     !<
1317    INTEGER(iwp) ::  ibase  !<
1318    INTEGER(iwp) ::  ic     !<
1319    INTEGER(iwp) ::  id     !<
1320    INTEGER(iwp) ::  ie     !<
1321    INTEGER(iwp) ::  if     !<
1322    INTEGER(iwp) ::  igo    !<
1323    INTEGER(iwp) ::  iink   !<
1324    INTEGER(iwp) ::  ijk    !<
1325    INTEGER(iwp) ::  j      !<
1326    INTEGER(iwp) ::  ja     !<
1327    INTEGER(iwp) ::  jb     !<
1328    INTEGER(iwp) ::  jbase  !<
1329    INTEGER(iwp) ::  jc     !<
1330    INTEGER(iwp) ::  jd     !<
1331    INTEGER(iwp) ::  je     !<
1332    INTEGER(iwp) ::  jf     !<
1333    INTEGER(iwp) ::  jg     !<
1334    INTEGER(iwp) ::  jh     !<
1335    INTEGER(iwp) ::  jink   !<
1336    INTEGER(iwp) ::  jump   !<
1337    INTEGER(iwp) ::  k      !<
1338    INTEGER(iwp) ::  kb     !<
1339    INTEGER(iwp) ::  kc     !<
1340    INTEGER(iwp) ::  kd     !<
1341    INTEGER(iwp) ::  ke     !<
1342    INTEGER(iwp) ::  kf     !<
1343    INTEGER(iwp) ::  kstop  !<
1344    INTEGER(iwp) ::  l      !<
1345    INTEGER(iwp) ::  m      !<
1346
1347    REAL(wp) ::  a(*)      !<
1348    REAL(wp) ::  b(*)      !<
1349    REAL(wp) ::  c(*)      !<
1350    REAL(wp) ::  d(*)      !<
1351    REAL(wp) ::  trigs(*)  !<
1352
1353    REAL(wp) ::  c1      !<
1354    REAL(wp) ::  c2      !<
1355    REAL(wp) ::  c3      !<
1356    REAL(wp) ::  c4      !<
1357    REAL(wp) ::  c5      !<
1358    REAL(wp) ::  qqrt5   !<
1359    REAL(wp) ::  qrt5    !<
1360    REAL(wp) ::  s1      !<
1361    REAL(wp) ::  s2      !<
1362    REAL(wp) ::  s3      !<
1363    REAL(wp) ::  s4      !<
1364    REAL(wp) ::  s5      !<
1365    REAL(wp) ::  sin36   !<
1366    REAL(wp) ::  sin45   !<
1367    REAL(wp) ::  sin60   !<
1368    REAL(wp) ::  sin72   !<
1369    REAL(wp) ::  ssin36  !<
1370    REAL(wp) ::  ssin45  !<
1371    REAL(wp) ::  ssin60  !<
1372    REAL(wp) ::  ssin72  !<
1373
1374    REAL(wp) ::  a10(nfft)  !<
1375    REAL(wp) ::  a11(nfft)  !<
1376    REAL(wp) ::  a20(nfft)  !<
1377    REAL(wp) ::  a21(nfft)  !<
1378    REAL(wp) ::  b10(nfft)  !<
1379    REAL(wp) ::  b11(nfft)  !<
1380    REAL(wp) ::  b20(nfft)  !<
1381    REAL(wp) ::  b21(nfft)  !<
1382
1383
1384    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
1385          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
1386
1387
1388    ierr = 0
1389
1390    m = n / ifac
1391    iink = la * inc1
1392    jink = la * inc2
1393    jump = ( ifac - 1 ) * jink
1394    kstop = ( n - ifac ) / ( 2 * ifac )
1395
1396    IF ( lot > nfft )  THEN
1397       ierr = 1
1398       RETURN
1399    ENDIF
1400    ibase = 0
1401    jbase = 0
1402    igo = ifac - 1
1403    IF ( igo == 7 )  igo = 6
1404    IF ( igo < 1  .OR.  igo > 6 )  THEN
1405       ierr = 2
1406       RETURN
1407    ENDIF
1408
1409
1410    SELECT CASE ( igo )
1411!
1412!--    Coding for factor 2
1413       CASE ( 1 )
1414
1415          ia = 1
1416          ib = ia + ( 2 * m - la ) * inc1
1417          ja = 1
1418          jb = ja + jink
1419
1420          IF ( la /= m )  THEN
1421
1422             DO  l = 1, la
1423                i = ibase
1424                j = jbase
1425                !DIR$ IVDEP
1426                !CDIR NODEP
1427                !OCL NOVREC
1428                DO  ijk = 1, lot
1429                   c(ja+j) = a(ia+i) + a(ib+i)
1430                   c(jb+j) = a(ia+i) - a(ib+i)
1431                   i = i + inc3
1432                   j = j + inc4
1433                ENDDO
1434                ibase = ibase + inc1
1435                jbase = jbase + inc2
1436             ENDDO
1437
1438             ia = ia + iink
1439             iink  = 2 * iink
1440             ib = ib - iink
1441             ibase = 0
1442             jbase = jbase + jump
1443             jump = 2 * jump + jink
1444
1445             IF ( ia /= ib )  THEN
1446
1447                DO  k = la, kstop, la
1448                   kb = k + k
1449                   c1 = trigs(kb+1)
1450                   s1 = trigs(kb+2)
1451                   ibase = 0
1452                   DO  l = 1, la
1453                      i = ibase
1454                      j = jbase
1455                      !DIR$ IVDEP
1456                      !CDIR NODEP
1457                      !OCL NOVREC
1458                      DO  ijk = 1, lot
1459                         c(ja+j) = a(ia+i) + a(ib+i)
1460                         d(ja+j) = b(ia+i) - b(ib+i)
1461                         c(jb+j) = c1 * ( a(ia+i) - a(ib+i) ) - s1 * ( b(ia+i) + b(ib+i) )
1462                         d(jb+j) = s1 * ( a(ia+i) - a(ib+i) ) + c1 * ( b(ia+i) + b(ib+i) )
1463                         i = i + inc3
1464                         j = j + inc4
1465                      ENDDO
1466                      ibase = ibase + inc1
1467                      jbase = jbase + inc2
1468                   ENDDO
1469                   ia = ia + iink
1470                   ib = ib - iink
1471                   jbase = jbase + jump
1472                ENDDO
1473
1474                IF ( ia > ib )  RETURN
1475
1476             ENDIF
1477
1478             ibase = 0
1479             DO  l = 1, la
1480                i = ibase
1481                j = jbase
1482                !DIR$ IVDEP
1483                !CDIR NODEP
1484                !OCL NOVREC
1485                DO  ijk = 1, lot
1486                   c(ja+j) = a(ia+i)
1487                   c(jb+j) = - b(ia+i)
1488                   i = i + inc3
1489                   j = j + inc4
1490                ENDDO
1491                ibase = ibase + inc1
1492                jbase = jbase + inc2
1493             ENDDO
1494
1495          ELSE
1496
1497             DO  l = 1, la
1498                i = ibase
1499                j = jbase
1500                !DIR$ IVDEP
1501                !CDIR NODEP
1502                !OCL NOVREC
1503                DO  ijk = 1, lot
1504                   c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) )
1505                   c(jb+j) = 2.0_wp * ( a(ia+i) - a(ib+i) )
1506                   i = i + inc3
1507                   j = j + inc4
1508                ENDDO
1509                ibase = ibase + inc1
1510                jbase = jbase + inc2
1511             ENDDO
1512
1513          ENDIF
1514
1515!
1516!--    Coding for factor 3
1517       CASE ( 2 )
1518
1519          ia = 1
1520          ib = ia + ( 2 * m - la ) * inc1
1521          ic = ib
1522          ja = 1
1523          jb = ja + jink
1524          jc = jb + jink
1525
1526          IF ( la /= m )  THEN
1527
1528             DO  l = 1, la
1529                i = ibase
1530                j = jbase
1531                !DIR$ IVDEP
1532                !CDIR NODEP
1533                !OCL NOVREC
1534                DO  ijk = 1, lot
1535                   c(ja+j) = a(ia+i) + a(ib+i)
1536                   c(jb+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) - ( sin60 * ( b(ib+i) ) )
1537                   c(jc+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) + ( sin60 * ( b(ib+i) ) )
1538                   i = i + inc3
1539                   j = j + inc4
1540                ENDDO
1541                ibase = ibase + inc1
1542                jbase = jbase + inc2
1543             ENDDO
1544             ia = ia + iink
1545             iink = 2 * iink
1546             ib = ib + iink
1547             ic = ic - iink
1548             jbase = jbase + jump
1549             jump = 2 * jump + jink
1550
1551             IF ( ia /= ic )  THEN
1552
1553                DO  k = la, kstop, la
1554                   kb = k + k
1555                   kc = kb + kb
1556                   c1 = trigs(kb+1)
1557                   s1 = trigs(kb+2)
1558                   c2 = trigs(kc+1)
1559                   s2 = trigs(kc+2)
1560                   ibase = 0
1561                   DO  l = 1, la
1562                      i = ibase
1563                      j = jbase
1564                      !DIR$ IVDEP
1565                      !CDIR NODEP
1566                      !OCL NOVREC
1567                      DO  ijk = 1, lot
1568                         c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
1569                         d(ja+j) = b(ia+i) + ( b(ib+i) - b(ic+i) )
1570                         c(jb+j) = c1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60   &
1571                                   * ( b(ib+i) + b(ic+i) ) ) ) - s1 * ( ( b(ia+i) - 0.5_wp         &
1572                                   * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
1573                         d(jb+j) = s1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60   &
1574                                   * ( b(ib+i) + b(ic+i) ) ) ) + c1 * ( ( b(ia+i) - 0.5_wp         &
1575                                   * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
1576                         c(jc+j) = c2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60   &
1577                                   * ( b(ib+i) + b(ic+i) ) ) ) - s2 * ( ( b(ia+i) - 0.5_wp         &
1578                                   * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
1579                         d(jc+j) = s2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60   &
1580                                   * ( b(ib+i) + b(ic+i) ) ) ) + c2 * ( ( b(ia+i) - 0.5_wp         &
1581                                   * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
1582                         i = i + inc3
1583                         j = j + inc4
1584                      ENDDO
1585                      ibase = ibase + inc1
1586                      jbase = jbase + inc2
1587                   ENDDO
1588                   ia = ia + iink
1589                   ib = ib + iink
1590                   ic = ic - iink
1591                   jbase = jbase + jump
1592                ENDDO
1593
1594                IF ( ia > ic )  RETURN
1595
1596             ENDIF
1597
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)
1607                   c(jb+j) = ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) )
1608                   c(jc+j) = - ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) )
1609                   i = i + inc3
1610                   j = j + inc4
1611                ENDDO
1612                ibase = ibase + inc1
1613                jbase = jbase + inc2
1614             ENDDO
1615
1616          ELSE
1617
1618             ssin60 = 2.0_wp * sin60
1619             DO  l = 1, la
1620                i = ibase
1621                j = jbase
1622                !DIR$ IVDEP
1623                !CDIR NODEP
1624                !OCL NOVREC
1625                DO  ijk = 1, lot
1626                   c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) )
1627                   c(jb+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) - ( ssin60 * b(ib+i) )
1628                   c(jc+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) + ( ssin60 * b(ib+i) )
1629                   i = i + inc3
1630                   j = j + inc4
1631                ENDDO
1632                ibase = ibase + inc1
1633                jbase = jbase + inc2
1634             ENDDO
1635
1636          ENDIF
1637
1638!
1639!--    Coding for factor 4
1640       CASE ( 3 )
1641
1642          ia = 1
1643          ib = ia + ( 2 * m - la ) * inc1
1644          ic = ib + 2 * m * inc1
1645          id = ib
1646          ja = 1
1647          jb = ja + jink
1648          jc = jb + jink
1649          jd = jc + jink
1650
1651          IF ( la /= m )  THEN
1652
1653             DO  l = 1, la
1654                i = ibase
1655                j = jbase
1656                !DIR$ IVDEP
1657                !CDIR NODEP
1658                !OCL NOVREC
1659                DO  ijk = 1, lot
1660                   c(ja+j) = ( a(ia+i) + a(ic+i) ) + a(ib+i)
1661                   c(jb+j) = ( a(ia+i) - a(ic+i) ) - b(ib+i)
1662                   c(jc+j) = ( a(ia+i) + a(ic+i) ) - a(ib+i)
1663                   c(jd+j) = ( a(ia+i) - a(ic+i) ) + b(ib+i)
1664                   i = i + inc3
1665                   j = j + inc4
1666               ENDDO
1667                ibase = ibase + inc1
1668                jbase = jbase + inc2
1669             ENDDO
1670             ia    = ia + iink
1671             iink  = 2 * iink
1672             ib    = ib + iink
1673             ic    = ic - iink
1674             id    = id - iink
1675             jbase = jbase + jump
1676             jump  = 2 * jump + jink
1677
1678             IF ( ib /= ic )  THEN
1679
1680                DO  k = la, kstop, la
1681                   kb = k + k
1682                   kc = kb + kb
1683                   kd = kc + kb
1684                   c1 = trigs(kb+1)
1685                   s1 = trigs(kb+2)
1686                   c2 = trigs(kc+1)
1687                   s2 = trigs(kc+2)
1688                   c3 = trigs(kd+1)
1689                   s3 = trigs(kd+2)
1690                   ibase = 0
1691                   DO  l = 1, la
1692                      i = ibase
1693                      j = jbase
1694                      !DIR$ IVDEP
1695                      !CDIR NODEP
1696                      !OCL NOVREC
1697                      DO  ijk = 1, lot
1698                         c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) )
1699                         d(ja+j) = ( b(ia+i) - b(ic+i) ) + ( b(ib+i) - b(id+i) )
1700                         c(jc+j) = c2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) - s2     &
1701                                   * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) )
1702                         d(jc+j) = s2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) + c2     &
1703                                   * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) )
1704                         c(jb+j) = c1 * ( ( a(ia+i) - a(ic+i) ) - ( b(ib+i) + b(id+i) ) ) - s1     &
1705                                   * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) )
1706                         d(jb+j) = s1 * ( ( a(ia+i) - a(ic+i) ) - ( b (ib+i) + b(id+i) ) ) + c1    &
1707                                   * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) )
1708                         c(jd+j) = c3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) - s3     &
1709                                   * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) )
1710                         d(jd+j) = s3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) + c3     &
1711                                   * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) )
1712                         i = i + inc3
1713                         j = j + inc4
1714                      ENDDO
1715                      ibase = ibase + inc1
1716                      jbase = jbase + inc2
1717                   ENDDO
1718                   ia = ia + iink
1719                   ib = ib + iink
1720                   ic = ic - iink
1721                   id = id - iink
1722                   jbase = jbase + jump
1723                ENDDO
1724
1725                IF ( ib > ic )  RETURN
1726
1727             ENDIF
1728
1729             ibase = 0
1730             sin45 = SQRT( 0.5_wp )
1731             DO  l = 1, la
1732                i = ibase
1733                j = jbase
1734                !DIR$ IVDEP
1735                !CDIR NODEP
1736                !OCL NOVREC
1737                DO  ijk = 1, lot
1738                   c(ja+j) = a(ia+i) + a(ib+i)
1739                   c(jb+j) = sin45 * ( ( a(ia+i) - a(ib+i) ) - ( b(ia+i) + b(ib+i) ) )
1740                   c(jc+j) = b(ib+i) - b(ia+i)
1741                   c(jd+j) = - sin45 * ( ( a(ia+i) - a(ib+i) ) + ( b(ia+i) + b(ib+i) ) )
1742                   i = i + inc3
1743                   j = j + inc4
1744                ENDDO
1745                ibase = ibase + inc1
1746                jbase = jbase + inc2
1747             ENDDO
1748
1749          ELSE
1750
1751             DO  l = 1, la
1752                i = ibase
1753                j = jbase
1754                !DIR$ IVDEP
1755                !CDIR NODEP
1756                !OCL NOVREC
1757                DO  ijk = 1, lot
1758                   c(ja+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) + a(ib+i) )
1759                   c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) - b(ib+i) )
1760                   c(jc+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) - a(ib+i) )
1761                   c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) + b(ib+i) )
1762                   i = i + inc3
1763                   j = j + inc4
1764                ENDDO
1765                ibase = ibase + inc1
1766                jbase = jbase + inc2
1767             ENDDO
1768
1769          ENDIF
1770
1771!
1772!--    Coding for factor 5
1773       CASE ( 4 )
1774
1775          ia = 1
1776          ib = ia + ( 2 * m - la ) * inc1
1777          ic = ib + 2 * m * inc1
1778          id = ic
1779          ie = ib
1780          ja = 1
1781          jb = ja + jink
1782          jc = jb + jink
1783          jd = jc + jink
1784          je = jd + jink
1785
1786          IF ( la /= m )  THEN
1787
1788             DO  l = 1, la
1789                i = ibase
1790                j = jbase
1791                !DIR$ IVDEP
1792                !CDIR NODEP
1793                !OCL NOVREC
1794                DO  ijk = 1, lot
1795                   c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
1796                   c(jb+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 *              &
1797                             ( a(ib+i) - a(ic+i) ) ) - ( sin72 * b(ib+i) + sin36 * b(ic+i) )
1798                   c(jc+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 *              &
1799                             ( a(ib+i) - a(ic+i) ) ) - ( sin36 * b(ib+i) - sin72 * b(ic+i) )
1800                   c(jd+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 *              &
1801                             ( a(ib+i) - a(ic+i) ) ) + ( sin36 * b(ib+i) - sin72 * b(ic+i) )
1802                   c(je+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 *              &
1803                             ( a(ib+i) - a(ic+i) ) ) + ( sin72 * b(ib+i) + sin36 * b(ic+i) )
1804                   i = i + inc3
1805                   j = j + inc4
1806                ENDDO
1807                ibase = ibase + inc1
1808                jbase = jbase + inc2
1809             ENDDO
1810             ia    = ia + iink
1811             iink  = 2 * iink
1812             ib    = ib + iink
1813             ic    = ic + iink
1814             id    = id - iink
1815             ie    = ie - iink
1816             jbase = jbase + jump
1817             jump  = 2 * jump + jink
1818
1819             IF ( ib /= id )  THEN
1820
1821                DO  k = la, kstop, la
1822                   kb = k + k
1823                   kc = kb + kb
1824                   kd = kc + kb
1825                   ke = kd + kb
1826                   c1 = trigs(kb+1)
1827                   s1 = trigs(kb+2)
1828                   c2 = trigs(kc+1)
1829                   s2 = trigs(kc+2)
1830                   c3 = trigs(kd+1)
1831                   s3 = trigs(kd+2)
1832                   c4 = trigs(ke+1)
1833                   s4 = trigs(ke+2)
1834                   ibase = 0
1835                   DO  l = 1, la
1836                      i = ibase
1837                      j = jbase
1838                      !DIR$ IVDEP
1839                      !CDIR NODEP
1840                      !OCL NOVREC
1841                      DO  ijk = 1, lot
1842                         a10(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i)      &
1843                                    + a(id+i) ) ) ) + qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i)   &
1844                                    + a(id+i) ) )
1845                         a20(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i)      &
1846                                    + a(id+i) ) ) ) - qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i)   &
1847                                    + a(id+i) ) )
1848                         b10(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i)      &
1849                                    - b(id+i) ) ) ) + qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i)   &
1850                                    - b(id+i) ) )
1851                         b20(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i)      &
1852                                    - b(id+i) ) ) ) - qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i)   &
1853                                    - b(id+i) ) )
1854                         a11(ijk) = sin72 * ( b(ib+i) + b(ie+i) ) + sin36 * ( b(ic+i) + b(id+i) )
1855                         a21(ijk) = sin36 * ( b(ib+i) + b(ie+i) ) - sin72 * ( b(ic+i) + b(id+i) )
1856                         b11(ijk) = sin72 * ( a(ib+i) - a(ie+i) ) + sin36 * ( a(ic+i) - a(id+i) )
1857                         b21(ijk) = sin36 * ( a(ib+i) - a(ie+i) ) - sin72 * ( a(ic+i) - a(id+i) )
1858                         c(ja+j)  = a(ia+i) + ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i) + a(id+i) ) )
1859                         d(ja+j)  = b(ia+i) + ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(id+i) ) )
1860                         c(jb+j)  = c1 * ( a10(ijk) - a11(ijk) ) - s1 * ( b10(ijk) + b11(ijk) )
1861                         d(jb+j)  = s1 * ( a10(ijk) - a11(ijk) ) + c1 * ( b10(ijk) + b11(ijk) )
1862                         c(je+j)  = c4 * ( a10(ijk) + a11(ijk) ) - s4 * ( b10(ijk) - b11(ijk) )
1863                         d(je+j)  = s4 * ( a10(ijk) + a11(ijk) ) + c4 * ( b10(ijk) - b11(ijk) )
1864                         c(jc+j)  = c2 * ( a20(ijk) - a21(ijk) ) - s2 * ( b20(ijk) + b21(ijk) )
1865                         d(jc+j)  = s2 * ( a20(ijk) - a21(ijk) ) + c2 * ( b20(ijk) + b21(ijk) )
1866                         c(jd+j)  = c3 * ( a20(ijk) + a21(ijk) ) - s3 * ( b20(ijk) - b21(ijk) )
1867                         d(jd+j)  = s3 * ( a20(ijk) + a21(ijk) ) + c3 * ( b20(ijk) - b21(ijk) )
1868                         i = i + inc3
1869                         j = j + inc4
1870                      ENDDO
1871                      ibase = ibase + inc1
1872                      jbase = jbase + inc2
1873                   ENDDO
1874                   ia = ia + iink
1875                   ib = ib + iink
1876                   ic = ic + iink
1877                   id = id - iink
1878                   ie = ie - iink
1879                   jbase = jbase + jump
1880                ENDDO
1881
1882                IF ( ib > id )  RETURN
1883
1884             ENDIF
1885
1886             ibase = 0
1887             DO  l = 1, la
1888                i = ibase
1889                j = jbase
1890                !DIR$ IVDEP
1891                !CDIR NODEP
1892                !OCL NOVREC
1893                DO  ijk = 1, lot
1894                   c(ja+j) = ( a(ia+i) + a(ib+i) ) + a(ic+i)
1895                   c(jb+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) )    &
1896                             - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) )
1897                   c(je+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) )  &
1898                             - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) )
1899                   c(jc+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) )    &
1900                             - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) )
1901                   c(jd+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) )  &
1902                             - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) )
1903                   i = i + inc3
1904                   j = j + inc4
1905                ENDDO
1906                ibase = ibase + inc1
1907                jbase = jbase + inc2
1908             ENDDO
1909
1910          ELSE
1911
1912             qqrt5  = 2.0_wp * qrt5
1913             ssin36 = 2.0_wp * sin36
1914             ssin72 = 2.0_wp * sin72
1915             DO  l = 1, la
1916                i = ibase
1917                j = jbase
1918                !DIR$ IVDEP
1919                !CDIR NODEP
1920                !OCL NOVREC
1921                DO  ijk = 1, lot
1922                   c(ja+j) = 2.0_wp * ( a(ia+i) + ( a(ib+i) + a(ic+i) ) )
1923                   c(jb+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qqrt5      &
1924                             * ( a(ib+i) - a(ic+i) ) ) - ( ssin72 * b(ib+i) + ssin36 * b(ic+i) )
1925                   c(jc+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qqrt5      &
1926                             * ( a(ib+i) - a(ic+i) ) ) - ( ssin36 * b(ib+i) - ssin72 * b(ic+i) )
1927                   c(jd+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qqrt5      &
1928                             * ( a(ib+i) - a(ic+i) ) ) + ( ssin36 * b(ib+i) - ssin72 * b(ic+i) )
1929                   c(je+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qqrt5      &
1930                             * ( a(ib+i) - a(ic+i) ) ) + (ssin72 * b(ib+i) + ssin36 * b(ic+i) )
1931                   i = i + inc3
1932                   j = j + inc4
1933                ENDDO
1934                ibase = ibase + inc1
1935                jbase = jbase + inc2
1936             ENDDO
1937
1938          ENDIF
1939
1940!
1941!--    Coding for factor 6
1942       CASE ( 5 )
1943
1944          ia = 1
1945          ib = ia + ( 2 * m - la ) * inc1
1946          ic = ib + 2 * m * inc1
1947          id = ic + 2 * m * inc1
1948          ie = ic
1949          if = ib
1950          ja = 1
1951          jb = ja + jink
1952          jc = jb + jink
1953          jd = jc + jink
1954          je = jd + jink
1955          jf = je + jink
1956
1957          IF ( la /= m )  THEN
1958
1959             DO  l = 1, la
1960                i = ibase
1961                j = jbase
1962                !DIR$ IVDEP
1963                !CDIR NODEP
1964                !OCL NOVREC
1965                DO  ijk = 1, lot
1966                   c(ja+j) = ( a (ia+i) + a(id+i) ) +  ( a(ib+i) + a(ic+i) )
1967                   c(jd+j) = ( a (ia+i) - a(id+i) ) -  ( a(ib+i) - a(ic+i) )
1968                   c(jb+j) = ( ( a(ia+i) - a(id+i) ) + 0.5_wp * ( a(ib+i) - a(ic+i) ) )            &
1969                             - ( sin60 * ( b(ib+i) + b(ic+i) ) )
1970                   c(jf+j) = ( ( a(ia+i) - a(id+i) ) + 0.5_wp * ( a(ib+i) - a(ic+i) ) )            &
1971                             + ( sin60 * b(ib+i) + b(ic+i) )
1972                   c(jc+j) = ( ( a(ia+i) + a(id+i) ) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )            &
1973                             - ( sin60 * b(ib+i) - b(ic+i) )
1974                   c(je+j) = ( ( a(ia+i) + a(id+i) ) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )            &
1975                             + ( sin60 * b(ib+i) - b(ic+i) )
1976                   i = i + inc3
1977                   j = j + inc4
1978                ENDDO
1979                ibase = ibase + inc1
1980                jbase = jbase + inc2
1981             ENDDO
1982             ia    = ia + iink
1983             iink  = 2 * iink
1984             ib    = ib + iink
1985             ic    = ic + iink
1986             id    = id - iink
1987             ie    = ie - iink
1988             if    = if - iink
1989             jbase = jbase + jump
1990             jump  = 2 * jump + jink
1991
1992             IF ( ic /= id )  THEN
1993
1994                DO  k = la, kstop, la
1995                   kb = k + k
1996                   kc = kb + kb
1997                   kd = kc + kb
1998                   ke = kd + kb
1999                   kf = ke + kb
2000                   c1 = trigs(kb+1)
2001                   s1 = trigs(kb+2)
2002                   c2 = trigs(kc+1)
2003                   s2 = trigs(kc+2)
2004                   c3 = trigs(kd+1)
2005                   s3 = trigs(kd+2)
2006                   c4 = trigs(ke+1)
2007                   s4 = trigs(ke+2)
2008                   c5 = trigs(kf+1)
2009                   s5 = trigs(kf+2)
2010                   ibase = 0
2011                   DO  l = 1, la
2012                      i = ibase
2013                      j = jbase
2014                      !DIR$ IVDEP
2015                      !CDIR NODEP
2016                      !OCL NOVREC
2017                      DO  ijk = 1, lot
2018                         a11(ijk) = ( a(ie+i) + a(ib+i) ) + ( a(ic+i) + a(if+i) )
2019                         a20(ijk) = ( a(ia+i) + a(id+i) ) - 0.5_wp * a11(ijk)
2020                         a21(ijk) = sin60 * ( ( a(ie+i) + a(ib+i) ) - ( a(ic+i) + a(if+i) ) )
2021                         b11(ijk) = ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(if+i) )
2022                         b20(ijk) = ( b(ia+i) - b(id+i) ) - 0.5_wp * b11(ijk)
2023                         b21(ijk) = sin60 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i) - b(if+i) ) )
2024
2025                         c(ja+j) = ( a(ia+i) + a(id+i) ) + a11(ijk)
2026                         d(ja+j) = ( b(ia+i) - b(id+i) ) + b11(ijk)
2027                         c(jc+j) = c2 * ( a20(ijk) - b21(ijk) ) - s2 * ( b20(ijk) + a21(ijk) )
2028                         d(jc+j) = s2 * ( a20(ijk) - b21(ijk) ) + c2 * ( b20(ijk) + a21(ijk) )
2029                         c(je+j) = c4 * ( a20(ijk) + b21(ijk) ) - s4 * ( b20(ijk) - a21(ijk) )
2030                         d(je+j) = s4 * ( a20(ijk) + b21(ijk) ) + c4 * ( b20(ijk) - a21(ijk) )
2031
2032                         a11(ijk) = ( a(ie+i) - a(ib+i) ) + ( a(ic+i) - a(if+i) )
2033                         b11(ijk) = ( b(ie+i) + b(ib+i) ) - ( b(ic+i) + b(if+i) )
2034                         a20(ijk) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11(ijk)
2035                         a21(ijk) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
2036                         b20(ijk) = ( b(ia+i) + b(id+i) ) + 0.5_wp * b11(ijk)
2037                         b21(ijk) = sin60 * ( ( b(ie+i) + b(ib+i) ) + ( b(ic+i) + b(if+i) ) )
2038
2039                         c(jd+j) = c3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) - s3                   &
2040                                      * ( (b(ia+i) + b(id+i) ) - b11(ijk) )
2041                         d(jd+j) = s3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) + c3                   &
2042                                      * ( (b(ia+i) + b(id+i) ) - b11(ijk) )
2043                         c(jb+j) = c1 * ( a20(ijk) - b21(ijk) ) - s1 * ( b20(ijk) - a21(ijk) )
2044                         d(jb+j) = s1 * ( a20(ijk) - b21(ijk) ) + c1 * ( b20(ijk) - a21(ijk) )
2045                         c(jf+j) = c5 * ( a20(ijk) + b21(ijk) ) - s5 * ( b20(ijk) + a21(ijk) )
2046                         d(jf+j) = s5 * ( a20(ijk) + b21(ijk) ) + c5 * ( b20(ijk) + a21(ijk) )
2047
2048                         i = i + inc3
2049                         j = j + inc4
2050                      ENDDO
2051                      ibase = ibase + inc1
2052                      jbase = jbase + inc2
2053                   ENDDO
2054                   ia = ia + iink
2055                   ib = ib + iink
2056                   ic = ic + iink
2057                   id = id - iink
2058                   ie = ie - iink
2059                   if = if - iink
2060                   jbase = jbase + jump
2061                ENDDO
2062
2063                IF ( ic > id )  RETURN
2064
2065             ENDIF
2066
2067             ibase = 0
2068             DO  l = 1, la
2069                i = ibase
2070                j = jbase
2071                !DIR$ IVDEP
2072                !CDIR NODEP
2073                !OCL NOVREC
2074                DO  ijk = 1, lot
2075                   c(ja+j) = a(ib+i) + ( a(ia+i) + a(ic+i) )
2076                   c(jd+j) = b(ib+i) - ( b(ia+i) + b(ic+i) )
2077                   c(jb+j) = ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i) + b(ic+i) )  &
2078                             + b(ib+i) )
2079                   c(jf+j) = - ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i)            &
2080                             + b(ic+i) ) + b(ib+i) )
2081                   c(jc+j) = sin60 * ( b(ic+i) - b(ia+i) ) + ( 0.5_wp * ( a(ia+i) + a(ic+i) )      &
2082                             - a(ib+i) )
2083                   c(je+j) = sin60 * ( b(ic+i) - b(ia+i) ) - ( 0.5_wp * ( a(ia+i) + a(ic+i) )      &
2084                             - a(ib+i) )
2085                   i = i + inc3
2086                   j = j + inc4
2087                ENDDO
2088                ibase = ibase + inc1
2089                jbase = jbase + inc2
2090             ENDDO
2091
2092          ELSE
2093
2094             ssin60 = 2.0_wp * sin60
2095             DO  l = 1, la
2096                i = ibase
2097                j = jbase
2098                !DIR$ IVDEP
2099                !CDIR NODEP
2100                !OCL NOVREC
2101                DO  ijk = 1, lot
2102                   c(ja+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) ) + ( 2.0_wp * ( a(ib+i) + a(ic+i) ) )
2103                   c(jd+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) ) - ( 2.0_wp * ( a(ib+i) - a(ic+i) ) )
2104                   c(jb+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) + ( a(ib+i) - a(ic+i) ) )            &
2105                             - ( ssin60 * ( b(ib+i) + b(ic+i) ) )
2106                   c(jf+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) + ( a(ib+i) - a(ic+i) ) )            &
2107                             + ( ssin60 * ( b(ib+i) + b(ic+i) ) )
2108                   c(jc+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) - ( a(ib+i) + a(ic+i) ) )            &
2109                             - ( ssin60 * ( b(ib+i) - b(ic+i) ) )
2110                   c(je+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) - ( a(ib+i) + a(ic+i) ) )            &
2111                             + ( ssin60 * ( b(ib+i) - b(ic+i) ) )
2112                   i = i + inc3
2113                   j = j + inc4
2114                ENDDO
2115                ibase = ibase + inc1
2116                jbase = jbase + inc2
2117             ENDDO
2118
2119          ENDIF
2120
2121!
2122!--    Coding for factor 8
2123       CASE ( 6 )
2124
2125          IF ( la /= m )  THEN
2126             ierr = 3
2127             RETURN
2128          ENDIF
2129
2130          ia = 1
2131          ib = ia + la * inc1
2132          ic = ib + 2 * la * inc1
2133          id = ic + 2 * la * inc1
2134          ie = id + 2 * la * inc1
2135          ja = 1
2136          jb = ja + jink
2137          jc = jb + jink
2138          jd = jc + jink
2139          je = jd + jink
2140          jf = je + jink
2141          jg = jf + jink
2142          jh = jg + jink
2143          ssin45 = SQRT( 2.0_wp )
2144
2145          DO  l = 1, la
2146             i = ibase
2147             j = jbase
2148             !DIR$ IVDEP
2149             !CDIR NODEP
2150             !OCL NOVREC
2151             DO  ijk = 1, lot
2152                c(ja+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) )
2153                c(je+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) )
2154                c(jc+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) - ( b(ib+i) - b(id+i) ) )
2155                c(jg+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) + ( b(ib+i) - b(id+i) ) )
2156                c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) - b(ic+i) ) + ssin45 * ( ( a(ib+i)      &
2157                          - a(id+i) ) - ( b(ib+i) + b(id+i) ) )
2158                c(jf+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) - b(ic+i) ) - ssin45 * ( ( a(ib+i)      &
2159                          - a(id+i) ) - ( b(ib+i) + b(id+i) ) )
2160                c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) + b(ic+i) ) - ssin45 * ( ( a(ib+i)      &
2161                          - a(id+i) ) + ( b(ib+i) + b(id+i) ) )
2162                c(jh+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) + b(ic+i) ) + ssin45 * ( ( a(ib+i)      &
2163                          - a(id+i) ) + ( b(ib+i) + b(id+i) ) )
2164                i = i + inc3
2165                j = j + inc4
2166             ENDDO
2167             ibase = ibase + inc1
2168             jbase = jbase + inc2
2169          ENDDO
2170
2171    END SELECT
2172
2173 END SUBROUTINE rpassm
2174
2175!--------------------------------------------------------------------------------------------------!
2176! Description:
2177! ------------
2178!> Computes factors of n & trigonometric functins required by fft99 & fft991cy
2179!> Method: Dimension trigs(n),ifax(1),jfax(10),lfax(7)
2180!> Subroutine 'set99' - computes factors of n & trigonometric functins required by fft99 & fft991cy
2181!--------------------------------------------------------------------------------------------------!
2182 SUBROUTINE set99( trigs, ifax, n )
2183
2184
2185    USE control_parameters,                                                                        &
2186        ONLY:  message_string
2187
2188    USE kinds
2189
2190    IMPLICIT NONE
2191
2192    INTEGER(iwp) ::  ifax(*)   !<
2193
2194    INTEGER(iwp) ::  i         !<
2195    INTEGER(iwp) ::  ifac      !<
2196    INTEGER(iwp) ::  ixxx      !<
2197    INTEGER(iwp) ::  jfax(10)  !<
2198    INTEGER(iwp) ::  k         !<
2199    INTEGER(iwp) ::  l         !<
2200    INTEGER(iwp) ::  lfax(7)   !<
2201    INTEGER(iwp) ::  n         !<
2202    INTEGER(iwp) ::  nfax      !<
2203    INTEGER(iwp) ::  nhl       !<
2204    INTEGER(iwp) ::  nil       !<
2205    INTEGER(iwp) ::  nu        !<
2206
2207    REAL(wp) ::  angle         !<
2208    REAL(wp) ::  del           !<
2209
2210    REAL(wp) ::  trigs(*)      !<
2211
2212    DATA  lfax / 6, 8, 5, 4, 3, 2, 1 /
2213
2214
2215    ixxx = 1
2216
2217    del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND = wp )
2218    nil = 0
2219    nhl = ( n / 2 ) - 1
2220    DO  k = nil, nhl
2221       angle = REAL( k, KIND = wp ) * del
2222       trigs(2*k+1) = COS( angle )
2223       trigs(2*k+2) = SIN( angle )
2224    ENDDO
2225
2226!
2227!-- Find factors of n (8,6,5,4,3,2; only one 8 allowed).
2228!-- Look for sixes first, store factors in descending order.
2229    nu   = n
2230    ifac = 6
2231    k    = 0
2232    l    = 1
2233
2234    DO
2235
2236       IF ( MOD( nu, ifac ) == 0 )  THEN
2237          k = k + 1
2238          jfax(k) = ifac
2239          IF ( ifac == 8  .AND.  k /= 1 )  THEN
2240             jfax(1) = 8
2241             jfax(k) = 6
2242          ENDIF
2243          nu = nu/ifac
2244          IF ( nu == 1   )  EXIT
2245          IF ( ifac /= 8 )  CYCLE
2246       ENDIF
2247
2248       l = l + 1
2249       ifac = lfax(l)
2250       IF (ifac < 2 )  THEN
2251          message_string = 'number of gridpoints along x or/and y contain illegal  factors' //     &
2252                           '&only factors 2, 3, 5 are allowed'
2253          CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
2254          RETURN
2255       ENDIF
2256
2257    ENDDO
2258!
2259!-- Now reverse order of factors
2260    nfax = k
2261    ifax(1) = nfax
2262    DO  i = 1, nfax
2263       ifax(nfax+2-i) = jfax(i)
2264    ENDDO
2265    ifax(10) = n
2266
2267 END SUBROUTINE set99
2268
2269
2270 SUBROUTINE fft991cy_vec( a, work, trigs, ifax, n, isign )
2271
2272    USE kinds
2273
2274    IMPLICIT NONE
2275
2276
2277    INTEGER(iwp) ::  i       !<
2278    INTEGER(iwp) ::  ia      !<
2279    INTEGER(iwp) ::  ibase   !<
2280    INTEGER(iwp) ::  ierr    !<
2281    INTEGER(iwp) ::  ifac    !<
2282    INTEGER(iwp) ::  igo     !<
2283    INTEGER(iwp) ::  ii      !<
2284    INTEGER(iwp) ::  inc     !<
2285    INTEGER(iwp) ::  isign   !<
2286    INTEGER(iwp) ::  istart  !<
2287    INTEGER(iwp) ::  ix      !<
2288    INTEGER(iwp) ::  iz      !<
2289    INTEGER(iwp) ::  j       !<
2290    INTEGER(iwp) ::  jbase   !<
2291    INTEGER(iwp) ::  jump    !<
2292    INTEGER(iwp) ::  k       !<
2293    INTEGER(iwp) ::  la      !<
2294    INTEGER(iwp) ::  lot     !<
2295    INTEGER(iwp) ::  n       !<
2296    INTEGER(iwp) ::  nfax    !<
2297    INTEGER(iwp) ::  nvex    !<
2298    INTEGER(iwp) ::  nx      !<
2299
2300    INTEGER(iwp),DIMENSION(:),INTENT(IN) ::  ifax  !<
2301
2302    REAL(wp),DIMENSION(:,:) ::  a      !<
2303    REAL(wp),DIMENSION(:)   ::  trigs  !<
2304    REAL(wp),DIMENSION(:,:) ::  work   !<
2305
2306
2307    inc  = 1
2308    jump = n
2309    lot  = 1
2310
2311    IF ( ifax(10) /= n )  CALL set99( trigs, ifax, n )
2312    nfax = ifax(1)
2313    nx   = n + 1
2314    IF ( MOD( n, 2 ) == 1 )  nx = n
2315    nvex = 1
2316
2317    IF ( isign == 1 )  THEN
2318!
2319!--    Backward fft: spectral to gridpoint transform
2320
2321       istart = 1
2322       ia = istart
2323       i = istart
2324       a(:,i+inc) = 0.5_wp * a(:,i)
2325       IF ( MOD( n, 2 ) /= 1 )  THEN
2326          i = istart + n * inc
2327          a(:,i) = 0.5_wp * a(:,i)
2328       ENDIF
2329       ia = istart + inc
2330       la = 1
2331       igo = + 1
2332
2333       DO  k = 1, nfax
2334
2335          ifac = ifax(k+1)
2336          ierr = -1
2337          IF ( igo /= -1 )  THEN
2338             CALL rpassm_vec( a(:,ia:), a(:,ia+la*inc:), work, work(:,ifac*la+1:), trigs, inc, 1,  &
2339                              n, ifac, la, ierr )
2340          ELSE
2341             CALL rpassm_vec( work, work(:,la+1:), a(:,ia:), a(:,ia+ifac*la*inc:), trigs, 1, inc,  &
2342                              n, ifac, la, ierr )
2343          ENDIF
2344!
2345!--       Following messages shouldn't appear in PALM applications
2346          IF ( ierr /= 0 )  THEN
2347             SELECT CASE (ierr)
2348                CASE (:-1)
2349                   WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
2350                CASE (0)
2351                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
2352                CASE (1:)
2353                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
2354             END SELECT
2355             RETURN
2356          ENDIF
2357
2358          la  = ifac * la
2359          igo = -igo
2360          ia  = istart
2361
2362       ENDDO
2363!
2364!--    If necessary, copy results back to a
2365       IF ( MOD( nfax, 2 ) /= 0 )  THEN
2366          ibase = 1
2367          jbase = ia
2368          i = ibase
2369          j = jbase
2370          DO  ii = 1, n
2371             a(:,j) = work(:,i)
2372             i = i + 1
2373             j = j + inc
2374          ENDDO
2375       ENDIF
2376!
2377!--    Fill in zeros at end
2378       ix = istart + n * inc
2379       a(:,ix) = 0.0_wp
2380       a(:,ix+inc) = 0.0_wp
2381
2382    ELSEIF ( isign == -1 )  THEN
2383!
2384!--    Forward fft: gridpoint to spectral transform
2385       istart = 1
2386       ia  = istart
2387       la  = n
2388       igo = + 1
2389
2390       DO  k = 1, nfax
2391
2392          ifac = ifax(nfax+2-k)
2393          la = la / ifac
2394          ierr = -1
2395
2396          IF ( igo /= -1 )  THEN
2397             CALL qpassm_vec( a(:,ia:), a(:,ia+ifac*la*inc:), work, work(:,la+1:), trigs, inc, 1,  &
2398                              n, ifac, la, ierr )
2399          ELSE
2400             CALL qpassm_vec( work, work(:,ifac*la+1:), a(:,ia:), a(:,ia+la*inc:), trigs, 1, inc,  &
2401                              n, ifac, la, ierr )
2402          ENDIF
2403
2404!
2405!--       Following messages shouldn't appear in PALM applications
2406          IF ( ierr /= 0 )  THEN
2407             SELECT CASE (ierr)
2408                CASE (0)
2409                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
2410                CASE (1:)
2411                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
2412             END SELECT
2413             RETURN
2414          ENDIF
2415
2416          igo = -igo
2417          ia = istart + inc
2418
2419       ENDDO
2420!
2421!--    If necessary, copy results back to a
2422       IF ( MOD( nfax, 2 ) /= 0 )  THEN
2423          ibase = 1
2424          jbase = ia
2425          i = ibase
2426          j = jbase
2427          DO  ii = 1, n
2428             a(:,j) = work(:,i)
2429             i = i + 1
2430             j = j + inc
2431          ENDDO
2432       ENDIF
2433!
2434!--    Shift a(0) and fill in zero imag parts
2435       ix = istart
2436       a(:,ix) = a(:,ix+inc)
2437       a(:,ix+inc) = 0.0_wp
2438
2439       IF ( MOD( n, 2 ) /= 1 )  THEN
2440          iz = istart + ( n + 1 ) * inc
2441          a(:,iz) = 0.0_wp
2442       ENDIF
2443
2444    ENDIF
2445
2446 END SUBROUTINE fft991cy_vec
2447
2448!--------------------------------------------------------------------------------------------------!
2449! Description:
2450! ------------
2451!> Performs one pass through data as part of multiple real fft (fourier analysis) routine.
2452!>
2453!> Method:
2454!>
2455!> a       is first real input vector
2456!>         equivalence b(1) with a(ifac*la*inc1+1)
2457!> c       is first real output vector
2458!>         equivalence d(1) with c(la*inc2+1)
2459!> trigs   is a precalculated list of sines & cosines
2460!> inc1    is the addressing increment for a
2461!> inc2    is the addressing increment for c
2462!> inc3    is the increment between input vectors a
2463!> inc4    is the increment between output vectors c
2464!> lot     is the number of vectors
2465!> n       is the length of the vectors
2466!> ifac    is the current factor of n
2467!>         la = n/(product of factors used so far)
2468!> ierr    is an error indicator:
2469!>         0 - pass completed without error
2470!>         1 - lot greater than nfft
2471!>         2 - ifac not catered for
2472!>         3 - ifac only catered for if la=n/ifac
2473!--------------------------------------------------------------------------------------------------!
2474 SUBROUTINE qpassm_vec( a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr )
2475
2476    USE kinds
2477
2478    IMPLICIT NONE
2479
2480    INTEGER(iwp),INTENT(IN)  ::  ifac  !<
2481    INTEGER(iwp),INTENT(IN)  ::  inc1  !<
2482    INTEGER(iwp),INTENT(IN)  ::  inc2  !<
2483    INTEGER(iwp),INTENT(IN)  ::  la    !<
2484    INTEGER(iwp),INTENT(IN)  ::  n     !<
2485
2486    INTEGER(iwp),INTENT(OUT) ::  ierr  !<
2487
2488    INTEGER(iwp) ::  i      !<
2489    INTEGER(iwp) ::  ia     !<
2490    INTEGER(iwp) ::  ib     !<
2491    INTEGER(iwp) ::  ibase  !<
2492    INTEGER(iwp) ::  ic     !<
2493    INTEGER(iwp) ::  id     !<
2494    INTEGER(iwp) ::  ie     !<
2495    INTEGER(iwp) ::  if     !<
2496    INTEGER(iwp) ::  ig     !<
2497    INTEGER(iwp) ::  igo    !<
2498    INTEGER(iwp) ::  ih     !<
2499    INTEGER(iwp) ::  iink   !<
2500    INTEGER(iwp) ::  ijump  !<
2501    INTEGER(iwp) ::  ipl    !<  loop index parallel loop
2502    INTEGER(iwp) ::  j      !<
2503    INTEGER(iwp) ::  ja     !<
2504    INTEGER(iwp) ::  jb     !<
2505    INTEGER(iwp) ::  jbase  !<
2506    INTEGER(iwp) ::  jc     !<
2507    INTEGER(iwp) ::  jd     !<
2508    INTEGER(iwp) ::  je     !<
2509    INTEGER(iwp) ::  jf     !<
2510    INTEGER(iwp) ::  jink   !<
2511    INTEGER(iwp) ::  k      !<
2512    INTEGER(iwp) ::  kb     !<
2513    INTEGER(iwp) ::  kc     !<
2514    INTEGER(iwp) ::  kd     !<
2515    INTEGER(iwp) ::  ke     !<
2516    INTEGER(iwp) ::  kf     !<
2517    INTEGER(iwp) ::  kstop  !<
2518    INTEGER(iwp) ::  l      !<
2519    INTEGER(iwp) ::  m      !<
2520
2521!
2522!-- Arrays are dimensioned with n
2523    REAL(wp),DIMENSION(:,:) ::  a  !<
2524    REAL(wp),DIMENSION(:,:) ::  b  !<
2525    REAL(wp),DIMENSION(:,:) ::  c  !<
2526    REAL(wp),DIMENSION(:,:) ::  d  !<
2527
2528    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs  !<
2529
2530    REAL(wp) ::  a0      !<
2531    REAL(wp) ::  a1      !<
2532    REAL(wp) ::  a10     !<
2533    REAL(wp) ::  a11     !<
2534    REAL(wp) ::  a2      !<
2535    REAL(wp) ::  a20     !<
2536    REAL(wp) ::  a21     !<
2537    REAL(wp) ::  a3      !<
2538    REAL(wp) ::  a4      !<
2539    REAL(wp) ::  a5      !<
2540    REAL(wp) ::  a6      !<
2541    REAL(wp) ::  b0      !<
2542    REAL(wp) ::  b1      !<
2543    REAL(wp) ::  b10     !<
2544    REAL(wp) ::  b11     !<
2545    REAL(wp) ::  b2      !<
2546    REAL(wp) ::  b20     !<
2547    REAL(wp) ::  b21     !<
2548    REAL(wp) ::  b3      !<
2549    REAL(wp) ::  b4      !<
2550    REAL(wp) ::  b5      !<
2551    REAL(wp) ::  b6      !<
2552    REAL(wp) ::  c1      !<
2553    REAL(wp) ::  c2      !<
2554    REAL(wp) ::  c3      !<
2555    REAL(wp) ::  c4      !<
2556    REAL(wp) ::  c5      !<
2557    REAL(wp) ::  qrt5    !<
2558    REAL(wp) ::  s1      !<
2559    REAL(wp) ::  s2      !<
2560    REAL(wp) ::  s3      !<
2561    REAL(wp) ::  s4      !<
2562    REAL(wp) ::  s5      !<
2563    REAL(wp) ::  sin36   !<
2564    REAL(wp) ::  sin45   !<
2565    REAL(wp) ::  sin60   !<
2566    REAL(wp) ::  sin72   !<
2567    REAL(wp) ::  z       !<
2568    REAL(wp) ::  zqrt5   !<
2569    REAL(wp) ::  zsin36  !<
2570    REAL(wp) ::  zsin45  !<
2571    REAL(wp) ::  zsin60  !<
2572    REAL(wp) ::  zsin72  !<
2573
2574
2575
2576    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
2577          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
2578
2579
2580    ierr = 0
2581
2582    m = n / ifac
2583    iink  = la * inc1
2584    jink  = la * inc2
2585    ijump = ( ifac - 1 ) * iink
2586    kstop = ( n - ifac ) / ( 2 * ifac )
2587
2588    ibase = 0
2589    jbase = 0
2590    igo = ifac - 1
2591    IF ( igo == 7 )  igo = 6
2592
2593    IF (igo < 1 .OR. igo > 6 )  THEN
2594       ierr = 2
2595       RETURN
2596    ENDIF
2597
2598
2599    SELECT CASE ( igo )
2600
2601!
2602!--    Coding for factor 2
2603       CASE ( 1 )
2604          ia = 1
2605          ib = ia + iink
2606          ja = 1
2607          jb = ja + ( 2 * m - la ) * inc2
2608
2609          IF ( la /= m )  THEN
2610
2611             DO  l = 1, la
2612                i         = ibase
2613                j         = jbase
2614                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
2615                c(:,jb+j) = a(:,ia+i) - a(:,ib+i)
2616                ibase     = ibase + inc1
2617                jbase     = jbase + inc2
2618             ENDDO
2619
2620             ja    = ja + jink
2621             jink  = 2 * jink
2622             jb    = jb - jink
2623             ibase = ibase + ijump
2624             ijump = 2 * ijump + iink
2625
2626             IF ( ja /= jb )  THEN
2627
2628                DO  k = la, kstop, la
2629                   kb = k + k
2630                   c1 = trigs(kb+1)
2631                   s1 = trigs(kb+2)
2632                   jbase = 0
2633                   DO  l = 1, la
2634                      i = ibase
2635                      j = jbase
2636                      c(:,ja+j) = a(:,ia+i) + ( c1 * a(:,ib+i) + s1 * b(:,ib+i) )
2637                      c(:,jb+j) = a(:,ia+i) - ( c1 * a(:,ib+i) + s1 * b(:,ib+i) )
2638                      d(:,ja+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) + b(:,ia+i)
2639                      d(:,jb+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) - b(:,ia+i)
2640                      ibase = ibase + inc1
2641                      jbase = jbase + inc2
2642                   ENDDO
2643
2644                   ibase = ibase + ijump
2645                   ja = ja + jink
2646                   jb = jb - jink
2647                ENDDO
2648
2649                IF ( ja > jb )  RETURN
2650
2651             ENDIF
2652
2653             jbase = 0
2654             DO  l = 1, la
2655                i = ibase
2656                j = jbase
2657                c(:,ja+j) = a(:,ia+i)
2658                d(:,ja+j) = -a(:,ib+i)
2659                ibase = ibase + inc1
2660                jbase = jbase + inc2
2661             ENDDO
2662
2663          ELSE
2664
2665             z = 1.0_wp / REAL( n, KIND = wp )
2666             DO  l = 1, la
2667                i = ibase
2668                j = jbase
2669                c(:,ja+j) = z * ( a(:,ia+i) + a(:,ib+i) )
2670                c(:,jb+j) = z * ( a(:,ia+i) - a(:,ib+i) )
2671                ibase = ibase + inc1
2672                jbase = jbase + inc2
2673             ENDDO
2674
2675          ENDIF
2676
2677!
2678!--    Coding for factor 3
2679       CASE ( 2 )
2680
2681          ia = 1
2682          ib = ia + iink
2683          ic = ib + iink
2684          ja = 1
2685          jb = ja + ( 2 * m - la ) * inc2
2686          jc = jb
2687
2688          IF ( la /= m )  THEN
2689
2690             DO  l = 1, la
2691                i         = ibase
2692                j         = jbase
2693                c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
2694                c(:,jb+j) = a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) )
2695                d(:,jb+j) = sin60 * ( a(:,ic+i) - a(:,ib+i) )
2696                ibase     = ibase + inc1
2697                jbase     = jbase + inc2
2698             ENDDO
2699
2700             ja    = ja + jink
2701             jink  = 2 * jink
2702             jb    = jb + jink
2703             jc    = jc - jink
2704             ibase = ibase + ijump
2705             ijump = 2 * ijump + iink
2706
2707             IF ( ja /= jc )  THEN
2708
2709                DO  k = la, kstop, la
2710                   kb = k + k
2711                   kc = kb + kb
2712                   c1 = trigs(kb+1)
2713                   s1 = trigs(kb+2)
2714                   c2 = trigs(kc+1)
2715                   s2 = trigs(kc+2)
2716
2717                   jbase = 0
2718                   DO  l = 1, la
2719                      i = ibase
2720                      j = jbase
2721                      DO  ipl = 1, SIZE( a, 1 )
2722                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
2723                            + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
2724                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
2725                            + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
2726                         a2 = a(ipl,ia+i) - 0.5_wp * a1
2727                         b2 = b(ipl,ia+i) - 0.5_wp * b1
2728                         a3 = sin60 * ( ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                    &
2729                                      - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) )
2730                         b3 = sin60 * ( ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                    &
2731                                      - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) )
2732                         c(ipl,ja+j) = a(ipl,ia+i) + a1
2733                         d(ipl,ja+j) = b(ipl,ia+i) + b1
2734                         c(ipl,jb+j) = a2 + b3
2735                         d(ipl,jb+j) = b2 - a3
2736                         c(ipl,jc+j) = a2 - b3
2737                         d(ipl,jc+j) = - ( b2 + a3 )
2738                      ENDDO
2739                      ibase = ibase + inc1
2740                      jbase = jbase + inc2
2741                   ENDDO
2742                   ibase = ibase + ijump
2743                   ja    = ja + jink
2744                   jb    = jb + jink
2745                   jc    = jc - jink
2746                ENDDO
2747
2748                IF ( ja > jc )  RETURN
2749
2750             ENDIF
2751
2752             jbase = 0
2753             DO  l = 1, la
2754                i = ibase
2755                j = jbase
2756                c(:,ja+j) = a(:,ia+i) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) )
2757                d(:,ja+j) = - sin60 * ( a(:,ib+i) + a(:,ic+i) )
2758                c(:,jb+j) = a(:,ia+i) - ( a(:,ib+i) - a(:,ic+i) )
2759                ibase = ibase + inc1
2760                jbase = jbase + inc2
2761             ENDDO
2762
2763          ELSE
2764
2765             z = 1.0_wp / REAL( n, KIND = wp )
2766             zsin60 = z * sin60
2767             DO  l = 1, la
2768                i = ibase
2769                j = jbase
2770                c(:,ja+j) = z * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) )
2771                c(:,jb+j) = z * ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )
2772                d(:,jb+j) = zsin60 * ( a(:,ic+i) - a(:,ib+i) )
2773                ibase = ibase + inc1
2774                jbase = jbase + inc2
2775             ENDDO
2776
2777          ENDIF
2778
2779!
2780!--    Coding for factor 4
2781       CASE ( 3 )
2782
2783          ia = 1
2784          ib = ia + iink
2785          ic = ib + iink
2786          id = ic + iink
2787          ja = 1
2788          jb = ja + ( 2 * m - la ) * inc2
2789          jc = jb + 2 * m * inc2
2790          jd = jb
2791
2792          IF ( la /= m )  THEN
2793
2794             DO  l = 1, la
2795                i = ibase
2796                j = jbase
2797                c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) )
2798                c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) )
2799                c(:,jb+j) = a(:,ia+i) - a(:,ic+i)
2800                d(:,jb+j) = a(:,id+i) - a(:,ib+i)
2801                ibase = ibase + inc1
2802                jbase = jbase + inc2
2803             ENDDO
2804
2805             ja    = ja + jink
2806             jink  = 2 * jink
2807             jb    = jb + jink
2808             jc    = jc - jink
2809             jd    = jd - jink
2810             ibase = ibase + ijump
2811             ijump = 2 * ijump + iink
2812
2813             IF ( jb /= jc )  THEN
2814
2815                DO  k = la, kstop, la
2816                   kb = k + k
2817                   kc = kb + kb
2818                   kd = kc + kb
2819                   c1 = trigs(kb+1)
2820                   s1 = trigs(kb+2)
2821                   c2 = trigs(kc+1)
2822                   s2 = trigs(kc+2)
2823                   c3 = trigs(kd+1)
2824                   s3 = trigs(kd+2)
2825                   jbase = 0
2826                   DO  l = 1, la
2827                      i = ibase
2828                      j = jbase
2829                      DO  ipl = 1, SIZE( a, 1 )
2830                         a0 = a(ipl,ia+i) + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
2831                         a2 = a(ipl,ia+i) - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
2832                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
2833                            + ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) )
2834                         a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
2835                            - ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) )
2836                         b0 = b(ipl,ia+i) + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
2837                         b2 = b(ipl,ia+i) - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
2838                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
2839                            + ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) )
2840                         b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
2841                            - ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) )
2842                         c(ipl,ja+j) = a0 + a1
2843                         c(ipl,jc+j) = a0 - a1
2844                         d(ipl,ja+j) = b0 + b1
2845                         d(ipl,jc+j) = b1 - b0
2846                         c(ipl,jb+j) = a2 + b3
2847                         c(ipl,jd+j) = a2 - b3
2848                         d(ipl,jb+j) = b2 - a3
2849                         d(ipl,jd+j) = - ( b2 + a3 )
2850                      ENDDO
2851                      ibase = ibase + inc1
2852                      jbase = jbase + inc2
2853                   ENDDO
2854                   ibase = ibase + ijump
2855                   ja = ja + jink
2856                   jb = jb + jink
2857                   jc = jc - jink
2858                   jd = jd - jink
2859                ENDDO
2860
2861                IF ( jb > jc )  RETURN
2862
2863             ENDIF
2864
2865             sin45 = SQRT( 0.5_wp )
2866             jbase = 0
2867             DO  l = 1, la
2868                i = ibase
2869                j = jbase
2870                c(:,ja+j) =   a(:,ia+i) + sin45 * ( a(:,ib+i) - a(:,id+i) )
2871                c(:,jb+j) =   a(:,ia+i) - sin45 * ( a(:,ib+i) - a(:,id+i) )
2872                d(:,ja+j) = - a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) )
2873                d(:,jb+j) =   a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) )
2874                ibase = ibase + inc1
2875                jbase = jbase + inc2
2876             ENDDO
2877
2878          ELSE
2879
2880             z = 1.0_wp / REAL( n, KIND = wp )
2881             DO  l = 1, la
2882                i = ibase
2883                j = jbase
2884                c(:,ja+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) ) )
2885                c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )
2886                c(:,jb+j) = z * ( a(:,ia+i) - a(:,ic+i) )
2887                d(:,jb+j) = z * ( a(:,id+i) - a(:,ib+i) )
2888                ibase = ibase + inc1
2889                jbase = jbase + inc2
2890             ENDDO
2891
2892          ENDIF
2893
2894!
2895!--    Coding for factor 5
2896       CASE ( 4 )
2897
2898          ia = 1
2899          ib = ia + iink
2900          ic = ib + iink
2901          id = ic + iink
2902          ie = id + iink
2903          ja = 1
2904          jb = ja + ( 2 * m - la ) * inc2
2905          jc = jb + 2 * m * inc2
2906          jd = jc
2907          je = jb
2908
2909          IF ( la /= m )  THEN
2910
2911             DO  l = 1, la
2912                i = ibase
2913                j = jbase
2914                DO  ipl = 1, SIZE( a, 1 )
2915                   a1 = a(ipl,ib+i) + a(ipl,ie+i)
2916                   a3 = a(ipl,ib+i) - a(ipl,ie+i)
2917                   a2 = a(ipl,ic+i) + a(ipl,id+i)
2918                   a4 = a(ipl,ic+i) - a(ipl,id+i)
2919                   a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 )
2920                   a6 = qrt5 * ( a1 - a2 )
2921                   c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 )
2922                   c(ipl,jb+j) = a5 + a6
2923                   c(ipl,jc+j) = a5 - a6
2924                   d(ipl,jb+j) = - sin72 * a3 - sin36 * a4
2925                   d(ipl,jc+j) = - sin36 * a3 + sin72 * a4
2926                ENDDO
2927                ibase = ibase + inc1
2928                jbase = jbase + inc2
2929             ENDDO
2930
2931             ja    = ja + jink
2932             jink  = 2 * jink
2933             jb    = jb + jink
2934             jc    = jc + jink
2935             jd    = jd - jink
2936             je    = je - jink
2937             ibase = ibase + ijump
2938             ijump = 2 * ijump + iink
2939
2940             IF ( jb /= jd )  THEN
2941
2942                DO  k = la, kstop, la
2943                   kb = k + k
2944                   kc = kb + kb
2945                   kd = kc + kb
2946                   ke = kd + kb
2947                   c1 = trigs(kb+1)
2948                   s1 = trigs(kb+2)
2949                   c2 = trigs(kc+1)
2950                   s2 = trigs(kc+2)
2951                   c3 = trigs(kd+1)
2952                   s3 = trigs(kd+2)
2953                   c4 = trigs(ke+1)
2954                   s4 = trigs(ke+2)
2955                   jbase = 0
2956                   DO  l = 1, la
2957                      i = ibase
2958                      j = jbase
2959                      DO  ipl = 1, SIZE(a,1)
2960                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) + ( c4 * a(ipl,ie+i) + s4    &
2961                              * b(ipl,ie+i) )
2962                         a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) - ( c4 * a(ipl,ie+i) + s4    &
2963                              * b(ipl,ie+i) )
2964                         a2 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) + ( c3 * a(ipl,id+i) + s3    &
2965                              * b(ipl,id+i) )
2966                         a4 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) - ( c3 * a(ipl,id+i) + s3    &
2967                              * b(ipl,id+i) )
2968
2969                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) + ( c4 * b(ipl,ie+i) - s4    &
2970                              * a(ipl,ie+i) )
2971                         b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) - ( c4 * b(ipl,ie+i) - s4    &
2972                              * a(ipl,ie+i) )
2973                         b2 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) + ( c3 * b(ipl,id+i) - s3    &
2974                              * a(ipl,id+i) )
2975                         b4 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) - ( c3 * b(ipl,id+i) - s3    &
2976                              * a(ipl,id+i) )
2977
2978                         a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 )
2979                         a6 = qrt5 * ( a1 - a2 )
2980                         b5 = b(ipl,ia+i) - 0.25_wp * ( b1 + b2 )
2981                         b6 = qrt5 * ( b1 - b2)
2982
2983                         a10 = a5 + a6
2984                         a20 = a5 - a6
2985                         b10 = b5 + b6
2986                         b20 = b5 - b6
2987                         a11 = sin72 * b3 + sin36 * b4
2988                         a21 = sin36 * b3 - sin72 * b4
2989
2990                         b11 = sin72 * a3 + sin36 * a4
2991                         b21 = sin36 * a3 - sin72 * a4
2992
2993                         c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 )
2994                         c(ipl,jb+j) = a10 + a11
2995                         c(ipl,je+j) = a10 - a11
2996                         c(ipl,jc+j) = a20 + a21
2997                         c(ipl,jd+j) = a20 - a21
2998
2999                         d(ipl,ja+j) = b(ipl,ia+i) + ( b1 + b2 )
3000                         d(ipl,jb+j) = b10 - b11
3001                         d(ipl,je+j) = - ( b10 + b11 )
3002                         d(ipl,jc+j) = b20 - b21
3003                         d(ipl,jd+j) = - ( b20 + b21 )
3004                      ENDDO
3005                      ibase = ibase + inc1
3006                      jbase = jbase + inc2
3007                   ENDDO
3008                   ibase = ibase + ijump
3009                   ja    = ja + jink
3010                   jb    = jb + jink
3011                   jc    = jc + jink
3012                   jd    = jd - jink
3013                   je    = je - jink
3014                ENDDO
3015
3016                IF ( jb > jd )  RETURN
3017
3018             ENDIF
3019
3020             jbase = 0
3021             DO  l = 1, la
3022                i = ibase
3023                j = jbase
3024                DO  ipl = 1, SIZE( a, 1 )
3025                   a1 = a(ipl,ib+i) + a(ipl,ie+i)
3026                   a3 = a(ipl,ib+i) - a(ipl,ie+i)
3027                   a2 = a(ipl,ic+i) + a(ipl,id+i)
3028                   a4 = a(ipl,ic+i) - a(ipl,id+i)
3029                   a5 = a(ipl,ia+i) + 0.25_wp * ( a3 - a4 )
3030                   a6 = qrt5 * ( a3 + a4 )
3031
3032                   c(ipl,ja+j) = a5 + a6
3033                   c(ipl,jb+j) = a5 - a6
3034                   c(ipl,jc+j) = a(ipl,ia+i) - ( a3 - a4 )
3035
3036                   d(ipl,ja+j) = - sin36 * a1 - sin72 * a2
3037                   d(ipl,jb+j) = - sin72 * a1 + sin36 * a2
3038                ENDDO
3039                ibase = ibase + inc1
3040                jbase = jbase + inc2
3041             ENDDO
3042
3043          ELSE
3044
3045             z = 1.0_wp / REAL( n, KIND = wp )
3046             zqrt5  = z * qrt5
3047             zsin36 = z * sin36
3048             zsin72 = z * sin72
3049             DO  l = 1, la
3050                i = ibase
3051                j = jbase
3052                DO  ipl = 1, SIZE( a, 1 )
3053                   a1 = a(ipl,ib+i) + a(ipl,ie+i)
3054                   a3 = a(ipl,ib+i) - a(ipl,ie+i)
3055                   a2 = a(ipl,ic+i) + a(ipl,id+i)
3056                   a4 = a(ipl,ic+i) - a(ipl,id+i)
3057                   a5 = z * ( a(ipl,ia+i) - 0.25_wp * ( a1 + a2 ) )
3058                   a6 = zqrt5 * ( a1 - a2 )
3059
3060                   c(ipl,ja+j) = z * ( a(ipl,ia+i) + ( a1 + a2 ) )
3061                   c(ipl,jb+j) = a5 + a6
3062                   c(ipl,jc+j) = a5 - a6
3063
3064                   d(ipl,jb+j) = - zsin72 * a3 - zsin36 * a4
3065                   d(ipl,jc+j) = - zsin36 * a3 + zsin72 * a4
3066                ENDDO
3067                ibase = ibase + inc1
3068                jbase = jbase + inc2
3069             ENDDO
3070
3071          ENDIF
3072
3073!
3074!--    Coding for factor 6
3075       CASE ( 5 )
3076
3077          ia = 1
3078          ib = ia + iink
3079          ic = ib + iink
3080          id = ic + iink
3081          ie = id + iink
3082          if = ie + iink
3083          ja = 1
3084          jb = ja + ( 2 * m - la ) * inc2
3085          jc = jb + 2 * m * inc2
3086          jd = jc + 2 * m * inc2
3087          je = jc
3088          jf = jb
3089
3090          IF ( la /= m )  THEN
3091
3092             DO  l = 1, la
3093                i = ibase
3094                j = jbase
3095                DO  ipl = 1, SIZE( a, 1 )
3096                   a11         = ( a(ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) )
3097                   c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11
3098                   c(ipl,jc+j) = ( a(ipl,ia+i) + a(ipl,id+i) - 0.5_wp * a11 )
3099                   d(ipl,jc+j) = sin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) ) -                         &
3100                                 ( a(ipl,ib+i) + a(ipl,ie+i) ) )
3101                   a11         = ( a(ipl,ic+i) - a(ipl,if+i) ) + ( a(ipl,ie+i) - a(ipl,ib+i) )
3102                   c(ipl,jb+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11
3103                   d(ipl,jb+j) = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) -                         &
3104                                 ( a(ipl,ic+i) - a(ipl,if+i) ) )
3105                   c(ipl,jd+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) + a11
3106                END DO
3107                ibase = ibase + inc1
3108                jbase = jbase + inc2
3109             ENDDO
3110             ja = ja + jink
3111             jink = 2 * jink
3112             jb = jb + jink
3113             jc = jc + jink
3114             jd = jd - jink
3115             je = je - jink
3116             jf = jf - jink
3117             ibase = ibase + ijump
3118             ijump = 2 * ijump + iink
3119
3120             IF ( jc /= jd )  THEN
3121
3122                DO  k = la, kstop, la
3123                   kb = k + k
3124                   kc = kb + kb
3125                   kd = kc + kb
3126                   ke = kd + kb
3127                   kf = ke + kb
3128                   c1 = trigs(kb+1)
3129                   s1 = trigs(kb+2)
3130                   c2 = trigs(kc+1)
3131                   s2 = trigs(kc+2)
3132                   c3 = trigs(kd+1)
3133                   s3 = trigs(kd+2)
3134                   c4 = trigs(ke+1)
3135                   s4 = trigs(ke+2)
3136                   c5 = trigs(kf+1)
3137                   s5 = trigs(kf+2)
3138                   jbase = 0
3139                   DO  l = 1, la
3140                      i = ibase
3141                      j = jbase
3142                      DO  ipl = 1, SIZE( a, 1 )
3143                         a1 = c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i)
3144                         b1 = c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i)
3145                         a2 = c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i)
3146                         b2 = c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i)
3147                         a3 = c3 * a(ipl,id+i) + s3 * b(ipl,id+i)
3148                         b3 = c3 * b(ipl,id+i) - s3 * a(ipl,id+i)
3149                         a4 = c4 * a(ipl,ie+i) + s4 * b(ipl,ie+i)
3150                         b4 = c4 * b(ipl,ie+i) - s4 * a(ipl,ie+i)
3151                         a5 = c5 * a(ipl,if+i) + s5 * b(ipl,if+i)
3152                         b5 = c5 * b(ipl,if+i) - s5 * a(ipl,if+i)
3153                         a11 = ( a2 + a5 ) + ( a1 + a4 )
3154                         a20 = ( a(ipl,ia+i) + a3 ) - 0.5_wp * a11
3155                         a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) )
3156                         b11 = ( b2 + b5 ) + ( b1 + b4 )
3157                         b20 = ( b(ipl,ia+i) + b3 ) - 0.5_wp * b11
3158                         b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) )
3159                         c(ipl,ja+j) = ( a(ipl,ia+i) + a3 ) + a11
3160                         d(ipl,ja+j) = ( b(ipl,ia+i) + b3 ) + b11
3161                         c(ipl,jc+j) = a20 - b21
3162                         d(ipl,jc+j) = a21 + b20
3163                         c(ipl,je+j) = a20 + b21
3164                         d(ipl,je+j) = a21 - b20
3165                         a11 = ( a2 - a5 ) + ( a4 - a1 )
3166                         a20 = ( a(ipl,ia+i) - a3 ) - 0.5_wp * a11
3167                         a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) )
3168                         b11 = ( b5 - b2 ) - ( b4 - b1 )
3169                         b20 = ( b3 - b(ipl,ia+i) ) - 0.5_wp * b11
3170                         b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) )
3171                         c(ipl,jb+j) = a20 - b21
3172                         d(ipl,jb+j) = a21 - b20
3173                         c(ipl,jd+j) = a11 + ( a(ipl,ia+i) - a3 )
3174                         d(ipl,jd+j) = b11 + ( b3 - b(ipl,ia+i) )
3175                         c(ipl,jf+j) = a20 + b21
3176                         d(ipl,jf+j) = a21 + b20
3177                      ENDDO
3178                      ibase = ibase + inc1
3179                      jbase = jbase + inc2
3180                   ENDDO
3181                   ibase = ibase + ijump
3182                   ja = ja + jink
3183                   jb = jb + jink
3184                   jc = jc + jink
3185                   jd = jd - jink
3186                   je = je - jink
3187                   jf = jf - jink
3188                ENDDO
3189
3190                IF ( jc > jd )  RETURN
3191
3192             ENDIF
3193
3194             jbase = 0
3195             DO  l = 1, la
3196                i = ibase
3197                j = jbase
3198                c(:,ja+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) + sin60             &
3199                          * ( a(:,ib+i) - a(:,if+i) )
3200                d(:,ja+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) - sin60           &
3201                            * ( a(:,ic+i) + a(:,ie+i) )
3202                c(:,jb+j) = a(:,ia+i) - ( a(:,ic+i) - a(:,ie+i) )
3203                d(:,jb+j) = a(:,id+i) - ( a(:,ib+i) + a(:,if+i) )
3204                c(:,jc+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) - sin60             &
3205                          * ( a(:,ib+i) - a(:,if+i) )
3206                d(:,jc+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) + sin60           &
3207                            * ( a(:,ic+i) + a(:,ie+i) )
3208                ibase = ibase + inc1
3209                jbase = jbase + inc2
3210             ENDDO
3211
3212          ELSE
3213
3214             z = 1.0_wp / REAL( n, KIND = wp )
3215             zsin60 = z * sin60
3216             DO  l = 1, la
3217                i = ibase
3218                j = jbase
3219                DO  ipl = 1, SIZE( a, 1 )
3220                   a11 = ( a (ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) )
3221                   c(ipl,ja+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) + a11 )
3222                   c(ipl,jc+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11 )
3223                   d(ipl,jc+j) = zsin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) )                          &
3224                                          - ( a(ipl,ib+i) + a(ipl,ie+i) ) )
3225                   a11 = ( a(ipl,ic+i) - a(ipl,if+i) )  + ( a(ipl,ie+i) - a(ipl,ib+i) )
3226                   c(ipl,jb+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11 )
3227                   d(ipl,jb+j) = zsin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) )                          &
3228                                          - ( a(ipl,ic+i) - a(ipl,if+i) ) )
3229                   c(ipl,jd+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )
3230                ENDDO
3231                ibase = ibase + inc1
3232                jbase = jbase + inc2
3233             ENDDO
3234
3235          ENDIF
3236
3237!
3238!--    Coding for factor 8
3239       CASE ( 6 )
3240
3241          IF ( la /= m )  THEN
3242             ierr = 3
3243             RETURN
3244          ENDIF
3245
3246          ia = 1
3247          ib = ia + iink
3248          ic = ib + iink
3249          id = ic + iink
3250          ie = id + iink
3251          if = ie + iink
3252          ig = if + iink
3253          ih = ig + iink
3254          ja = 1
3255          jb = ja + la * inc2
3256          jc = jb + 2 * m * inc2
3257          jd = jc + 2 * m * inc2
3258          je = jd + 2 * m * inc2
3259          z = 1.0_wp / REAL( n, KIND = wp )
3260          zsin45 = z * SQRT( 0.5_wp )
3261
3262          DO  l = 1, la
3263             i = ibase
3264             j = jbase
3265             c(:,ja+j) = z * ( ( ( a(:,ia+i) + a(:,ie+i) ) + ( a(:,ic+i) + a(:,ig+i) ) )           &
3266                         + ( ( a(:,id+i) + a(:,ih+i) ) + ( a(:,ib+i) + a(:,if+i) ) ) )
3267             c(:,je+j) = z * ( ( ( a(:,ia+i) + a(:,ie+i) ) + ( a(:,ic+i) + a(:,ig+i) ) )           &
3268                         - ( ( a(:,id+i) + a(:,ih+i) ) + ( a(:,ib+i) + a(:,if+i) ) ) )
3269             c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ie+i) ) - ( a(:,ic+i) + a(:,ig+i) ) )
3270             d(:,jc+j) = z * ( ( a(:,id+i) + a(:,ih+i) ) - ( a(:,ib+i) + a(:,if+i) ) )
3271             c(:,jb+j) = z * ( a(:,ia+i) - a(:,ie+i) ) + zsin45 * ( ( a(:,ih+i) - a(:,id+i) )      &
3272                         - ( a(:,if+i) - a(:,ib+i) ) )
3273             c(:,jd+j) = z * ( a(:,ia+i) - a(:,ie+i) ) - zsin45 * ( ( a(:,ih+i) - a(:,id+i) )      &
3274                         - ( a(:,if+i) - a(:,ib+i) ) )
3275             d(:,jb+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) )        &
3276                         + z * ( a(:,ig+i) - a(:,ic+i) )
3277             d(:,jd+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) )        &
3278                         - z * ( a(:,ig+i) - a(:,ic+i) )
3279             ibase = ibase + inc1
3280             jbase = jbase + inc2
3281          ENDDO
3282
3283    END SELECT
3284
3285 END SUBROUTINE qpassm_vec
3286
3287!--------------------------------------------------------------------------------------------------!
3288! Description:
3289! ------------
3290!> Same as qpassm, but for backward fft
3291!--------------------------------------------------------------------------------------------------!
3292 SUBROUTINE rpassm_vec(a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr )
3293
3294    USE kinds
3295
3296    IMPLICIT NONE
3297
3298
3299    INTEGER(iwp) ::  i      !<
3300    INTEGER(iwp) ::  ia     !<
3301    INTEGER(iwp) ::  ib     !<
3302    INTEGER(iwp) ::  ibase  !<
3303    INTEGER(iwp) ::  ic     !<
3304    INTEGER(iwp) ::  id     !<
3305    INTEGER(iwp) ::  ie     !<
3306    INTEGER(iwp) ::  ierr   !<
3307    INTEGER(iwp) ::  ifac   !<
3308    INTEGER(iwp) ::  if     !<
3309    INTEGER(iwp) ::  igo    !<
3310    INTEGER(iwp) ::  iink   !<
3311    INTEGER(iwp) ::  inc1   !<
3312    INTEGER(iwp) ::  inc2   !<
3313    INTEGER(iwp) ::  ipl    !<  loop index parallel loop
3314    INTEGER(iwp) ::  j      !<
3315    INTEGER(iwp) ::  ja     !<
3316    INTEGER(iwp) ::  jb     !<
3317    INTEGER(iwp) ::  jbase  !<
3318    INTEGER(iwp) ::  jc     !<
3319    INTEGER(iwp) ::  jd     !<
3320    INTEGER(iwp) ::  je     !<
3321    INTEGER(iwp) ::  jf     !<
3322    INTEGER(iwp) ::  jg     !<
3323    INTEGER(iwp) ::  jh     !<
3324    INTEGER(iwp) ::  jink   !<
3325    INTEGER(iwp) ::  jump   !<
3326    INTEGER(iwp) ::  k      !<
3327    INTEGER(iwp) ::  kb     !<
3328    INTEGER(iwp) ::  kc     !<
3329    INTEGER(iwp) ::  kd     !<
3330    INTEGER(iwp) ::  ke     !<
3331    INTEGER(iwp) ::  kf     !<
3332    INTEGER(iwp) ::  kstop  !<
3333    INTEGER(iwp) ::  l      !<
3334    INTEGER(iwp) ::  la     !<
3335    INTEGER(iwp) ::  m      !<
3336    INTEGER(iwp) ::  n      !<
3337!
3338!-- Arrays are dimensioned with n
3339    REAL(wp),DIMENSION(:,:) ::  a  !<
3340    REAL(wp),DIMENSION(:,:) ::  b  !<
3341    REAL(wp),DIMENSION(:,:) ::  c  !<
3342    REAL(wp),DIMENSION(:,:) ::  d  !<
3343    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs  !<
3344
3345    REAL(wp) ::  a10     !<
3346    REAL(wp) ::  a11     !<
3347    REAL(wp) ::  a20     !<
3348    REAL(wp) ::  a21     !<
3349    REAL(wp) ::  b10     !<
3350    REAL(wp) ::  b11     !<
3351    REAL(wp) ::  b20     !<
3352    REAL(wp) ::  b21     !<
3353    REAL(wp) ::  c1      !<
3354    REAL(wp) ::  c2      !<
3355    REAL(wp) ::  c3      !<
3356    REAL(wp) ::  c4      !<
3357    REAL(wp) ::  c5      !<
3358    REAL(wp) ::  qqrt5   !<
3359    REAL(wp) ::  qrt5    !<
3360    REAL(wp) ::  s1      !<
3361    REAL(wp) ::  s2      !<
3362    REAL(wp) ::  s3      !<
3363    REAL(wp) ::  s4      !<
3364    REAL(wp) ::  s5      !<
3365    REAL(wp) ::  sin36   !<
3366    REAL(wp) ::  sin45   !<
3367    REAL(wp) ::  sin60   !<
3368    REAL(wp) ::  sin72   !<
3369    REAL(wp) ::  ssin36  !<
3370    REAL(wp) ::  ssin45  !<
3371    REAL(wp) ::  ssin60  !<
3372    REAL(wp) ::  ssin72  !<
3373
3374
3375
3376    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
3377          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
3378
3379
3380    ierr = 0
3381
3382    m = n / ifac
3383    iink = la * inc1
3384    jink = la * inc2
3385    jump = ( ifac - 1 ) * jink
3386    kstop = ( n - ifac ) / ( 2 * ifac )
3387
3388    ibase = 0
3389    jbase = 0
3390    igo = ifac - 1
3391    IF ( igo == 7 )  igo = 6
3392    IF ( igo < 1  .OR.  igo > 6 )  THEN
3393       ierr = 2
3394       RETURN
3395    ENDIF
3396
3397
3398    SELECT CASE ( igo )
3399!
3400!--    Coding for factor 2
3401       CASE ( 1 )
3402
3403          ia = 1
3404          ib = ia + ( 2 * m - la ) * inc1
3405          ja = 1
3406          jb = ja + jink
3407
3408          IF ( la /= m )  THEN
3409
3410             DO  l = 1, la
3411                i = ibase
3412                j = jbase
3413                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
3414                c(:,jb+j) = a(:,ia+i) - a(:,ib+i)
3415                ibase = ibase + inc1
3416                jbase = jbase + inc2
3417             ENDDO
3418
3419             ia = ia + iink
3420             iink = 2 * iink
3421             ib = ib - iink
3422             ibase = 0
3423             jbase = jbase + jump
3424             jump = 2 * jump + jink
3425
3426             IF ( ia /= ib )  THEN
3427
3428                DO  k = la, kstop, la
3429                   kb = k + k
3430                   c1 = trigs(kb+1)
3431                   s1 = trigs(kb+2)
3432                   ibase = 0
3433                   DO  l = 1, la
3434                      i = ibase
3435                      j = jbase
3436                      c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
3437                      d(:,ja+j) = b(:,ia+i) - b(:,ib+i)
3438                      c(:,jb+j) = c1 * ( a(:,ia+i) - a(:,ib+i) ) - s1 * ( b(:,ia+i) + b(:,ib+i) )
3439                      d(:,jb+j) = s1 * ( a(:,ia+i) - a(:,ib+i) ) + c1 * ( b(:,ia+i) + b(:,ib+i) )
3440                      ibase = ibase + inc1
3441                      jbase = jbase + inc2
3442                   ENDDO
3443                   ia = ia + iink
3444                   ib = ib - iink
3445                   jbase = jbase + jump
3446                ENDDO
3447
3448                IF ( ia > ib )  RETURN
3449
3450             ENDIF
3451
3452             ibase = 0
3453             DO  l = 1, la
3454                i = ibase
3455                j = jbase
3456                c(:,ja+j) = a(:,ia+i)
3457                c(:,jb+j) = - b(:,ia+i)
3458                ibase = ibase + inc1
3459                jbase = jbase + inc2
3460             ENDDO
3461
3462          ELSE
3463
3464             DO  l = 1, la
3465                i = ibase
3466                j = jbase
3467                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) )
3468                c(:,jb+j) = 2.0_wp * ( a(:,ia+i) - a(:,ib+i) )
3469                ibase = ibase + inc1
3470                jbase = jbase + inc2
3471             ENDDO
3472
3473          ENDIF
3474
3475!
3476!--    Coding for factor 3
3477       CASE ( 2 )
3478
3479          ia = 1
3480          ib = ia + ( 2 * m - la ) * inc1
3481          ic = ib
3482          ja = 1
3483          jb = ja + jink
3484          jc = jb + jink
3485
3486          IF ( la /= m )  THEN
3487
3488             DO  l = 1, la
3489                i = ibase
3490                j = jbase
3491                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
3492                c(:,jb+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) - ( sin60 * ( b(:,ib+i) ) )
3493                c(:,jc+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) + ( sin60 * ( b(:,ib+i) ) )
3494                ibase = ibase + inc1
3495                jbase = jbase + inc2
3496             ENDDO
3497             ia = ia + iink
3498             iink = 2 * iink
3499             ib = ib + iink
3500             ic = ic - iink
3501             jbase = jbase + jump
3502             jump = 2 * jump + jink
3503
3504             IF ( ia /= ic )  THEN
3505
3506                DO  k = la, kstop, la
3507                   kb = k + k
3508                   kc = kb + kb
3509                   c1 = trigs(kb+1)
3510                   s1 = trigs(kb+2)
3511                   c2 = trigs(kc+1)
3512                   s2 = trigs(kc+2)
3513                   ibase = 0
3514                   DO  l = 1, la
3515                      i = ibase
3516                      j = jbase
3517                      c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
3518                      d(:,ja+j) = b(:,ia+i) + ( b(:,ib+i) - b(:,ic+i) )
3519                      c(:,jb+j) = c1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
3520                                  - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s1 * ( ( b(:,ia+i)   &
3521                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i)   &
3522                                  - a(:,ic+i) ) ) )
3523                      d(:,jb+j) = s1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
3524                                  - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c1 * ( ( b(:,ia+i)   &
3525                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i)   &
3526                                  - a(:,ic+i) ) ) )
3527                      c(:,jc+j) = c2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
3528                                  + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s2 * ( ( b(:,ia+i)   &
3529                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i)   &
3530                                  - a(:,ic+i) ) ) )
3531                      d(:,jc+j) = s2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
3532                                  + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c2 * ( ( b(:,ia+i)   &
3533                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i)   &
3534                                  - a(:,ic+i) ) ) )
3535                      ibase = ibase + inc1
3536                      jbase = jbase + inc2
3537                   ENDDO
3538                   ia = ia + iink
3539                   ib = ib + iink
3540                   ic = ic - iink
3541                   jbase = jbase + jump
3542                ENDDO
3543
3544                IF ( ia > ic )  RETURN
3545
3546             ENDIF
3547
3548             ibase = 0
3549             DO  l = 1, la
3550                i = ibase
3551                j = jbase
3552                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
3553                c(:,jb+j) = ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) )
3554                c(:,jc+j) = - ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) )
3555                ibase = ibase + inc1
3556                jbase = jbase + inc2
3557             ENDDO
3558
3559          ELSE
3560
3561             ssin60 = 2.0_wp * sin60
3562             DO  l = 1, la
3563                i = ibase
3564                j = jbase
3565                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) )
3566                c(:,jb+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) - ( ssin60 * b(:,ib+i) )
3567                c(:,jc+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) + ( ssin60 * b(:,ib+i) )
3568                ibase = ibase + inc1
3569                jbase = jbase + inc2
3570             ENDDO
3571
3572          ENDIF
3573
3574!
3575!--    Coding for factor 4
3576       CASE ( 3 )
3577
3578          ia = 1
3579          ib = ia + ( 2 * m - la ) * inc1
3580          ic = ib + 2 * m * inc1
3581          id = ib
3582          ja = 1
3583          jb = ja + jink
3584          jc = jb + jink
3585          jd = jc + jink
3586
3587          IF ( la /= m )  THEN
3588
3589             DO  l = 1, la
3590                i = ibase
3591                j = jbase
3592                c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i)
3593                c(:,jb+j) = ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i)
3594                c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i)
3595                c(:,jd+j) = ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i)
3596                ibase = ibase + inc1
3597                jbase = jbase + inc2
3598             ENDDO
3599             ia = ia + iink
3600             iink = 2 * iink
3601             ib = ib + iink
3602             ic = ic - iink
3603             id = id - iink
3604             jbase = jbase + jump
3605             jump = 2 * jump + jink
3606
3607             IF ( ib /= ic )  THEN
3608
3609                DO  k = la, kstop, la
3610                   kb = k + k
3611                   kc = kb + kb
3612                   kd = kc + kb
3613                   c1 = trigs(kb+1)
3614                   s1 = trigs(kb+2)
3615                   c2 = trigs(kc+1)
3616                   s2 = trigs(kc+2)
3617                   c3 = trigs(kd+1)
3618                   s3 = trigs(kd+2)
3619                   ibase = 0
3620                   DO  l = 1, la
3621                      i = ibase
3622                      j = jbase
3623                      c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) )
3624                      d(:,ja+j) = ( b(:,ia+i) - b(:,ic+i) ) + ( b(:,ib+i) - b(:,id+i) )
3625                      c(:,jc+j) = c2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )   &
3626                                  - s2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) )
3627                      d(:,jc+j) = s2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )   &
3628                                  + c2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) )
3629                      c(:,jb+j) = c1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) )   &
3630                                  - s1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) )
3631                      d(:,jb+j) = s1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) )   &
3632                                  + c1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) )
3633                      c(:,jd+j) = c3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) )   &
3634                                  - s3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) )
3635                      d(:,jd+j) = s3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) )   &
3636                                  + c3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) )
3637                      ibase = ibase + inc1
3638                      jbase = jbase + inc2
3639                   ENDDO
3640                   ia = ia + iink
3641                   ib = ib + iink
3642                   ic = ic - iink
3643                   id = id - iink
3644                   jbase = jbase + jump
3645                ENDDO
3646
3647                IF ( ib > ic )  RETURN
3648
3649             ENDIF
3650
3651             ibase = 0
3652             sin45 = SQRT( 0.5_wp )
3653             DO  l = 1, la
3654                i = ibase
3655                j = jbase
3656                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
3657                c(:,jb+j) = sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) - ( b(:,ia+i) + b(:,ib+i) ) )
3658                c(:,jc+j) = b(:,ib+i) - b(:,ia+i)
3659                c(:,jd+j) = - sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) + ( b(:,ia+i) + b(:,ib+i) ) )
3660                ibase = ibase + inc1
3661                jbase = jbase + inc2
3662             ENDDO
3663
3664          ELSE
3665
3666             DO  l = 1, la
3667                i = ibase
3668                j = jbase
3669                c(:,ja+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i) )
3670                c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i) )
3671                c(:,jc+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i) )
3672                c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i) )
3673                ibase = ibase + inc1
3674                jbase = jbase + inc2
3675             ENDDO
3676
3677          ENDIF
3678
3679!
3680!--    Coding for factor 5
3681       CASE ( 4 )
3682
3683          ia = 1
3684          ib = ia + ( 2 * m - la ) * inc1
3685          ic = ib + 2 * m * inc1
3686          id = ic
3687          ie = ib
3688          ja = 1
3689          jb = ja + jink
3690          jc = jb + jink
3691          jd = jc + jink
3692          je = jd + jink
3693
3694          IF ( la /= m )  THEN
3695
3696             DO  l = 1, la
3697                i = ibase
3698                j = jbase
3699                c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
3700                c(:,jb+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5           &
3701                          * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) )
3702                c(:,jc+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5           &
3703                          * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) )
3704                c(:,jd+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5           &
3705                          * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) )
3706                c(:,je+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5           &
3707                          * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) )
3708                ibase = ibase + inc1
3709                jbase = jbase + inc2
3710             ENDDO
3711             ia = ia + iink
3712             iink = 2 * iink
3713             ib = ib + iink
3714             ic = ic + iink
3715             id = id - iink
3716             ie = ie - iink
3717             jbase = jbase + jump
3718             jump = 2 * jump + jink
3719
3720             IF ( ib /= id )  THEN
3721
3722                DO  k = la, kstop, la
3723                   kb = k + k
3724                   kc = kb + kb
3725                   kd = kc + kb
3726                   ke = kd + kb
3727                   c1 = trigs(kb+1)
3728                   s1 = trigs(kb+2)
3729                   c2 = trigs(kc+1)
3730                   s2 = trigs(kc+2)
3731                   c3 = trigs(kd+1)
3732                   s3 = trigs(kd+2)
3733                   c4 = trigs(ke+1)
3734                   s4 = trigs(ke+2)
3735                   ibase = 0
3736                   DO  l = 1, la
3737                      i = ibase
3738                      j = jbase
3739!DIR$ IVDEP
3740                      DO  ipl = 1, SIZE( a, 1 )
3741                         a10 = ( a (ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) )          &
3742                             + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) + qrt5 * ( ( a(ipl,ib+i)          &
3743                             + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) )
3744                         a20 = ( a(ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) )           &
3745                             + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) - qrt5 * ( ( a(ipl,ib+i)          &
3746                             + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) )
3747                         b10 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) )           &
3748                             + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) + qrt5 * ( ( b(ipl,ib+i)          &
3749                             - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) )
3750                         b20 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) )           &
3751                             + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) - qrt5 * ( ( b(ipl,ib+i)          &
3752                             - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) )
3753                         a11 = sin72 * ( b(ipl,ib+i) + b(ipl,ie+i) ) + sin36                       &
3754                               * ( b(ipl,ic+i) + b(ipl,id+i) )
3755                         a21 = sin36 * ( b(ipl,ib+i) + b(ipl,ie+i) ) - sin72                       &
3756                               * ( b(ipl,ic+i) + b(ipl,id+i) )
3757                         b11 = sin72 * ( a(ipl,ib+i) - a(ipl,ie+i) ) + sin36                       &
3758                               * ( a(ipl,ic+i) - a(ipl,id+i) )
3759                         b21 = sin36 * ( a(ipl,ib+i) - a(ipl,ie+i) ) - sin72                       &
3760                               * ( a(ipl,ic+i) - a(ipl,id+i) )
3761                         c(ipl,ja+j) = a(ipl,ia+i) + ( ( a(ipl,ib+i) + a(ipl,ie+i) )               &
3762                                       + ( a(ipl,ic+i) + a(ipl,id+i) ) )
3763                         d(ipl,ja+j) = b(ipl,ia+i) + ( ( b(ipl,ib+i) - b(ipl,ie+i) )               &
3764                                       + ( b(ipl,ic+i) - b(ipl,id+i) ) )
3765                         c(ipl,jb+j) = c1 * ( a10 - a11 ) - s1 * ( b10 + b11 )
3766                         d(ipl,jb+j) = s1 * ( a10 - a11 ) + c1 * ( b10 + b11 )
3767                         c(ipl,je+j) = c4 * ( a10 + a11 ) - s4 * ( b10 - b11 )
3768                         d(ipl,je+j) = s4 * ( a10 + a11 ) + c4 * ( b10 - b11 )
3769                         c(ipl,jc+j) = c2 * ( a20 - a21 ) - s2 * ( b20 + b21 )
3770                         d(ipl,jc+j) = s2 * ( a20 - a21 ) + c2 * ( b20 + b21 )
3771                         c(ipl,jd+j) = c3 * ( a20 + a21 ) - s3 * ( b20 - b21 )
3772                         d(ipl,jd+j) = s3 * ( a20 + a21 ) + c3 * ( b20 - b21 )
3773                      ENDDO
3774                      ibase = ibase + inc1
3775                      jbase = jbase + inc2
3776                   ENDDO
3777                   ia = ia + iink
3778                   ib = ib + iink
3779                   ic = ic + iink
3780                   id = id - iink
3781                   ie = ie - iink
3782                   jbase = jbase + jump
3783                ENDDO
3784
3785                IF ( ib > id )  RETURN
3786
3787             ENDIF
3788
3789             ibase = 0
3790             DO  l = 1, la
3791                i = ibase
3792                j = jbase
3793                c(:,ja+j) = ( a(:,ia+i) + a(:,ib+i) ) + a(:,ic+i)
3794                c(:,jb+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i)           &
3795                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) )
3796                c(:,je+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i)         &
3797                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) )
3798                c(:,jc+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i)           &
3799                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) )
3800                c(:,jd+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i)         &
3801                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) )
3802                ibase = ibase + inc1
3803                jbase = jbase + inc2
3804             ENDDO
3805
3806          ELSE
3807
3808             qqrt5  = 2.0_wp * qrt5
3809             ssin36 = 2.0_wp * sin36
3810             ssin72 = 2.0_wp * sin72
3811             DO  l = 1, la
3812                i = ibase
3813                j = jbase
3814                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) )
3815                c(:,jb+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
3816                            + qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) )  - ( ssin72 * b(:,ib+i)          &
3817                            + ssin36 * b(:,ic+i) )
3818                c(:,jc+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
3819                            - qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) - ( ssin36 * b(:,ib+i)           &
3820                            - ssin72 * b(:,ic+i) )
3821                c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
3822                            - qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( ssin36 * b(:,ib+i)           &
3823                            - ssin72 * b(:,ic+i) )
3824                c(:,je+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
3825                            + qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( ssin72 * b(:,ib+i)           &
3826                            + ssin36 * b(:,ic+i) )
3827                ibase = ibase + inc1
3828                jbase = jbase + inc2
3829             ENDDO
3830
3831          ENDIF
3832
3833!
3834!--    Coding for factor 6
3835       CASE ( 5 )
3836
3837          ia = 1
3838          ib = ia + ( 2 * m - la ) * inc1
3839          ic = ib + 2 * m * inc1
3840          id = ic + 2 * m * inc1
3841          ie = ic
3842          if = ib
3843          ja = 1
3844          jb = ja + jink
3845          jc = jb + jink
3846          jd = jc + jink
3847          je = jd + jink
3848          jf = je + jink
3849
3850          IF ( la /= m )  THEN
3851
3852             DO  l = 1, la
3853                i = ibase
3854                j = jbase
3855                c(:,ja+j) = ( a(:,ia+i) + a(:,id+i) ) +  ( a(:,ib+i) + a(:,ic+i) )
3856                c(:,jd+j) = ( a(:,ia+i) - a(:,id+i) ) -  ( a(:,ib+i) - a(:,ic+i) )
3857                c(:,jb+j) = ( ( a(:,ia+i) - a(:,id+i) ) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) )     &
3858                          - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) )
3859                c(:,jf+j) = ( ( a(:,ia+i) - a(:,id+i) ) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) )     &
3860                          + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) )
3861                c(:,jc+j) = ( ( a(:,ia+i) + a(:,id+i) ) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )     &
3862                          - ( sin60 * ( b(:,ib+i) - b(:,ic+i) ) )
3863                c(:,je+j) = ( ( a(:,ia+i) + a(:,id+i) ) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )     &
3864                          + ( sin60 * ( b(:,ib+i) - b(:,ic+i) ) )
3865                ibase = ibase + inc1
3866                jbase = jbase + inc2
3867             ENDDO
3868             ia = ia + iink
3869             iink = 2 * iink
3870             ib = ib + iink
3871             ic = ic + iink
3872             id = id - iink
3873             ie = ie - iink
3874             if = if - iink
3875             jbase = jbase + jump
3876             jump = 2 * jump + jink
3877
3878             IF ( ic /= id )  THEN
3879
3880                DO  k = la, kstop, la
3881                   kb = k + k
3882                   kc = kb + kb
3883                   kd = kc + kb
3884                   ke = kd + kb
3885                   kf = ke + kb
3886                   c1 = trigs(kb+1)
3887                   s1 = trigs(kb+2)
3888                   c2 = trigs(kc+1)
3889                   s2 = trigs(kc+2)
3890                   c3 = trigs(kd+1)
3891                   s3 = trigs(kd+2)
3892                   c4 = trigs(ke+1)
3893                   s4 = trigs(ke+2)
3894                   c5 = trigs(kf+1)
3895                   s5 = trigs(kf+2)
3896                   ibase = 0
3897                   DO  l = 1, la
3898                      i = ibase
3899                      j = jbase
3900                      DO  ipl = 1, SIZE( a, 1)
3901                         a11 = ( a(ipl,ie+i) + a(ipl,ib+i) ) + ( a(ipl,ic+i) + a(ipl,if+i) )
3902                         a20 = ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11
3903                         a21 = sin60 * ( ( a(ipl,ie+i) + a(ipl,ib+i) )                             &
3904                               - ( a(ipl,ic+i) + a(ipl,if+i) ) )
3905                         b11 = ( b(ipl,ib+i) - b(ipl,ie+i) ) + ( b(ipl,ic+i) - b(ipl,if+i) )
3906                         b20 = ( b(ipl,ia+i) - b(ipl,id+i) ) - 0.5_wp * b11
3907                         b21 = sin60 * ( ( b(ipl,ib+i) - b(ipl,ie+i) )                             &
3908                               - ( b(ipl,ic+i) - b(ipl,if+i) ) )
3909
3910                         c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11
3911                         d(ipl,ja+j) = ( b(ipl,ia+i) - b(ipl,id+i) ) + b11
3912                         c(ipl,jc+j) = c2 * ( a20 - b21 ) - s2 * ( b20 + a21 )
3913                         d(ipl,jc+j) = s2 * ( a20 - b21 ) + c2 * ( b20 + a21 )
3914                         c(ipl,je+j) = c4 * ( a20 + b21 ) - s4 * ( b20 - a21 )
3915                         d(ipl,je+j) = s4 * ( a20 + b21 ) + c4 * ( b20 - a21 )
3916
3917                         a11 = ( a(ipl,ie+i) - a(ipl,ib+i) ) + ( a(ipl,ic+i) - a(ipl,if+i) )
3918                         b11 = ( b(ipl,ie+i) + b(ipl,ib+i) ) - ( b(ipl,ic+i) + b(ipl,if+i) )
3919                         a20 = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11
3920                         a21 = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) - ( a(ipl,ic+i)             &
3921                               - a(ipl,if+i) ) )
3922                         b20 = ( b(ipl,ia+i) + b(ipl,id+i) ) + 0.5_wp * b11
3923                         b21 = sin60 * ( ( b(ipl,ie+i) + b(ipl,ib+i) ) + ( b(ipl,ic+i)             &
3924                               + b(ipl,if+i) ) )
3925
3926                         c(ipl,jd+j) = c3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )                &
3927                                     - s3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 )
3928                         d(ipl,jd+j) = s3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )                &
3929                                     + c3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 )
3930                         c(ipl,jb+j) = c1 * ( a20 - b21 ) - s1 * ( b20 - a21 )
3931                         d(ipl,jb+j) = s1 * ( a20 - b21 ) + c1 * ( b20 - a21 )
3932                         c(ipl,jf+j) = c5 * ( a20 + b21 ) - s5 * ( b20 + a21 )
3933                         d(ipl,jf+j) = s5 * ( a20 + b21 ) + c5 * ( b20 + a21 )
3934
3935                      ENDDO
3936                      ibase = ibase + inc1
3937                      jbase = jbase + inc2
3938                   ENDDO
3939                   ia = ia + iink
3940                   ib = ib + iink
3941                   ic = ic + iink
3942                   id = id - iink
3943                   ie = ie - iink
3944                   if = if - iink
3945                   jbase = jbase + jump
3946                ENDDO
3947
3948                IF ( ic > id )  RETURN
3949
3950             ENDIF
3951
3952             ibase = 0
3953             DO  l = 1, la
3954                i = ibase
3955                j = jbase
3956                DO  ipl = 1, SIZE( a ,1 )
3957                   c(ipl,ja+j) = a(ipl,ib+i) + ( a(ipl,ia+i) + a(ipl,ic+i) )
3958                   c(ipl,jd+j) = b(ipl,ib+i) - ( b(ipl,ia+i) + b(ipl,ic+i) )
3959                   c(ipl,jb+j) = ( sin60 * ( a(ipl,ia+i) - a(ipl,ic+i) ) )                         &
3960                               - ( 0.5_wp * ( b(ipl,ia+i) + b(ipl,ic+i) ) + b(ipl,ib+i) )
3961                   c(ipl,jf+j) = - ( sin60 * ( a(ipl,ia+i) - a(ipl,ic+i) ) )                       &
3962                                 - ( 0.5_wp * ( b(ipl,ia+i) + b(ipl,ic+i) ) + b(ipl,ib+i) )
3963                   c(ipl,jc+j) = sin60 * ( b(ipl,ic+i) - b(ipl,ia+i) )                              &
3964                                 + ( 0.5_wp * ( a(ipl,ia+i) + a(ipl,ic+i) ) - a(ipl,ib+i) )
3965                   c(ipl,je+j) = sin60 * ( b(ipl,ic+i) - b(ipl,ia+i) )                              &
3966                                 - ( 0.5_wp * ( a(ipl,ia+i) + a(ipl,ic+i) ) - a(ipl,ib+i) )
3967                ENDDO
3968                ibase = ibase + inc1
3969                jbase = jbase + inc2
3970             ENDDO
3971
3972          ELSE
3973
3974             ssin60 = 2.0_wp * sin60
3975             DO  l = 1, la
3976                i = ibase
3977                j = jbase
3978                c(:,ja+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) )                                 &
3979                          + ( 2.0_wp * ( a(:,ib+i) + a(:,ic+i) ) )
3980                c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) )                                 &
3981                          - ( 2.0_wp * ( a(:,ib+i) - a(:,ic+i) ) )
3982                c(:,jb+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) + ( a(:,ib+i) - a(:,ic+i) ) )     &
3983                          - ( ssin60 * ( b(:,ib+i) + b(:,ic+i) ) )
3984                c(:,jf+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) + ( a(:,ib+i) - a(:,ic+i) ) )     &
3985                          + ( ssin60 * ( b(:,ib+i) + b(:,ic+i) ) )
3986                c(:,jc+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) - ( a(:,ib+i) + a(:,ic+i) ) )     &
3987                          - ( ssin60 * ( b(:,ib+i) - b(:,ic+i) ) )
3988                c(:,je+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) - ( a(:,ib+i) + a(:,ic+i) ) )     &
3989                          + ( ssin60 * ( b(:,ib+i) - b(:,ic+i) ) )
3990                ibase = ibase + inc1
3991                jbase = jbase + inc2
3992             ENDDO
3993
3994          ENDIF
3995
3996!
3997!--    Coding for factor 8
3998       CASE ( 6 )
3999
4000          IF ( la /= m )  THEN
4001             ierr = 3
4002             RETURN
4003          ENDIF
4004
4005          ia = 1
4006          ib = ia + la * inc1
4007          ic = ib + 2 * la * inc1
4008          id = ic + 2 * la * inc1
4009          ie = id + 2 * la * inc1
4010          ja = 1
4011          jb = ja + jink
4012          jc = jb + jink
4013          jd = jc + jink
4014          je = jd + jink
4015          jf = je + jink
4016          jg = jf + jink
4017          jh = jg + jink
4018          ssin45 = SQRT( 2.0_wp )
4019
4020          DO  l = 1, la
4021             i = ibase
4022             j = jbase
4023             c(:,ja+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) )                      &
4024                                                  + ( a(:,ib+i) + a(:,id+i) ) )
4025             c(:,je+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) )                      &
4026                                                - ( a(:,ib+i) + a(:,id+i) ) )
4027             c(:,jc+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) )                      &
4028                                                  - ( b(:,ib+i) - b(:,id+i) ) )
4029             c(:,jg+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) )                      &
4030                                                  + ( b(:,ib+i) - b(:,id+i) ) )
4031             c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) - b(:,ic+i) ) + ssin45               &
4032                                * ( (a(:,ib+i) - a(:,id+i) ) - ( b(:,ib+i) + b(:,id+i) ) )
4033             c(:,jf+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) - b(:,ic+i) ) - ssin45               &
4034                                * ( (a(:,ib+i) - a(:,id+i) ) - ( b(:,ib+i) + b(:,id+i) ) )
4035             c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) + b(:,ic+i) ) - ssin45               &
4036                                * ( (a(:,ib+i) - a(:,id+i) ) + ( b(:,ib+i) + b(:,id+i) ) )
4037             c(:,jh+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) + b(:,ic+i) ) + ssin45               &
4038                                * ( (a(:,ib+i) - a(:,id+i) ) + ( b(:,ib+i) + b(:,id+i) ) )
4039             ibase = ibase + inc1
4040             jbase = jbase + inc2
4041          ENDDO
4042
4043    END SELECT
4044
4045 END SUBROUTINE rpassm_vec
4046
4047 END MODULE temperton_fft
Note: See TracBrowser for help on using the repository browser.