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

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

code vectorization for NEC Aurora: vectorized version of Temperton FFT, vectorization of Newtor iteration for calculating the Obukhov length

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