source: palm/trunk/SOURCE/temperton_fft_mod.f90 @ 4370

Last change on this file since 4370 was 4370, checked in by raasch, 5 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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