Home

Context Navigation

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

Last change on this file since 4096 was 3725, checked in by raasch, 6 years ago
modifications to avoid compiler warnings about unused variables, temperton-fft: GOTO statements replaced, file re-formatted corresponding to coding standards, ssh-calls for compilations on remote systems modified to avoid output of login messages on specific systems changed again (palmbuild, reverted as before r3549), error messages for failed restarts extended (palmrun)
Property svn:keywords set to `Id`
File size: 71.2 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 3725 2019-02-07 10:11:02Z motisi $
11	! GOTO statements replaced, file re-formatted corresponding to coding standards
12	!
13	! 3241 2018-09-12 15:02:00Z raasch
14	! unused variables removed
15	!
16	! 3049 2018-05-29 13:52:36Z Giersch
17	! Error messages revised
18	!
19	! 3045 2018-05-28 07:55:41Z Giersch
20	! Error message revised
21	!
22	! 2300 2017-06-29 13:31:14Z raasch
23	! NEC related CPP directives removed
24	!
25	! 1851 2016-04-08 13:32:50Z maronga
26	!
27	! 1850 2016-04-08 13:29:27Z maronga
28	! Module renamed
29	!
30	! 1682 2015-10-07 23:56:08Z knoop
31	! Code annotations made doxygen readable
32	!
33	! 1342 2014-03-26 17:04:47Z kanani
34	! REAL constants defined as wp-kind
35	!
36	! 1322 2014-03-20 16:38:49Z raasch
37	! REAL constants defined as wp-kind
38	!
39	! 1320 2014-03-20 08:40:49Z raasch
40	! ONLY-attribute added to USE-statements,
41	! kind-parameters added to all INTEGER and REAL declaration statements,
42	! kinds are defined in new module kinds,
43	! old module precision_kind is removed,
44	! revision history before 2012 removed,
45	! comment fields (!:) to be used for variable explanations added to
46	! all variable declaration statements
47	!
48	! Revision 1.1 2003/03/12 16:41:59 raasch
49	! Initial revision
50	!
51	!
52	! Description:
53	! ------------
54	!> Fast Fourier transformation developed by Clive Temperton, ECMWF.
55	!------------------------------------------------------------------------------!
56	MODULE temperton_fft
57
58	USE kinds
59
60	IMPLICIT NONE
61
62	PRIVATE
63
64	PUBLIC set99, fft991cy
65
66	INTEGER(iwp), PARAMETER :: nfft = 32 !< maximum length of calls to *fft
67	INTEGER(iwp), PARAMETER :: nout = 6 !< standard output stream
68
69	CONTAINS
70
71	!------------------------------------------------------------------------------!
72	! Description:
73	! ------------
74	!> Calls Fortran-versions of fft's.
75	!>
76	!> Method:
77	!>
78	!> Subroutine 'fft991cy' - multiple fast real periodic transform
79	!> supersedes previous routine 'fft991cy'.
80	!>
81	!> Real transform of length n performed by removing redundant
82	!> operations from complex transform of length n.
83	!>
84	!> a is the array containing input & output data.
85	!> work is an area of size (n+1)*min(lot,nfft).
86	!> trigs is a previously prepared list of trig function values.
87	!> ifax is a previously prepared list of factors of n.
88	!> inc is the increment within each data 'vector'
89	!> (e.g. inc=1 for consecutively stored data).
90	!> jump is the increment between the start of each data vector.
91	!> n is the length of the data vectors.
92	!> lot is the number of data vectors.
93	!> isign = +1 for transform from spectral to gridpoint
94	!> = -1 for transform from gridpoint to spectral.
95	!>
96	!> ordering of coefficients:
97	!> a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
98	!> where b(0)=b(n/2)=0; (n+2) locations required.
99	!>
100	!> ordering of data:
101	!> x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
102	!>
103	!> Vectorization is achieved on cray by doing the transforms
104	!> in parallel.
105	!>
106	!> n must be composed of factors 2,3 & 5 but does not have to be even.
107	!>
108	!> definition of transforms:
109	!>
110	!> isign=+1: x(j)=sum(k=0,.,n-1)(c(k)exp(2ijk*pi/n))
111	!> where c(k)=a(k)+ib(k) and c(n-k)=a(k)-ib(k)
112	!>
113	!> isign=-1: a(k)=(1/n)sum(j=0,.,n-1)(x(j)cos(2jk*pi/n))
114	!> b(k)=-(1/n)sum(j=0,.,n-1)(x(j)sin(2jk*pi/n))
115	!>
116	!> calls fortran-versions of fft's !!!
117	!> dimension a(n),work(n),trigs(n),ifax(1)
118	!------------------------------------------------------------------------------!
119	SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
120
121	USE kinds
122
123	IMPLICIT NONE
124
125	INTEGER(iwp) :: inc !<
126	INTEGER(iwp) :: isign !<
127	INTEGER(iwp) :: jump !<
128	INTEGER(iwp) :: lot !<
129	INTEGER(iwp) :: n !<
130
131	REAL(wp) :: a(*) !<
132	REAL(wp) :: trigs(*) !<
133	REAL(wp) :: work(*) !<
134	INTEGER(iwp) :: ifax(*) !<
135
136	INTEGER(iwp) :: i !<
137	INTEGER(iwp) :: ia !<
138	INTEGER(iwp) :: ibase !<
139	INTEGER(iwp) :: ierr !<
140	INTEGER(iwp) :: ifac !<
141	INTEGER(iwp) :: igo !<
142	INTEGER(iwp) :: ii !<
143	INTEGER(iwp) :: istart !<
144	INTEGER(iwp) :: ix !<
145	INTEGER(iwp) :: iz !<
146	INTEGER(iwp) :: j !<
147	INTEGER(iwp) :: jbase !<
148	INTEGER(iwp) :: jj !<
149	INTEGER(iwp) :: k !<
150	INTEGER(iwp) :: la !<
151	INTEGER(iwp) :: nb !<
152	INTEGER(iwp) :: nblox !<
153	INTEGER(iwp) :: nfax !<
154	INTEGER(iwp) :: nvex !<
155	INTEGER(iwp) :: nx !<
156
157
158	IF ( ifax(10) /= n ) CALL set99( trigs, ifax, n )
159	nfax = ifax(1)
160	nx = n + 1
161	IF ( MOD(n,2) == 1 ) nx = n
162	nblox = 1 + (lot-1) / nfft
163	nvex = lot - (nblox-1) * nfft
164
165
166	IF ( isign == 1 ) THEN
167	!
168	!-- Backward fft: spectral to gridpoint transform
169
170	istart = 1
171	DO nb = 1, nblox
172
173	ia = istart
174	i = istart
175	!DIR$ IVDEP
176	!CDIR NODEP
177	!OCL NOVREC
178	DO j = 1, nvex
179	a(i+inc) = 0.5_wp * a(i)
180	i = i + jump
181	ENDDO
182	IF ( MOD(n,2) /= 1 ) THEN
183	i = istart + n * inc
184	DO j = 1, nvex
185	a(i) = 0.5_wp * a(i)
186	i = i + jump
187	ENDDO
188	ENDIF
189	ia = istart + inc
190	la = 1
191	igo = + 1
192
193	DO k = 1, nfax
194
195	ifac = ifax(k+1)
196	ierr = -1
197	IF ( igo /= -1 ) THEN
198	CALL rpassm( a(ia), a(ia+lainc), work(1), work(ifacla+1), trigs, inc, 1, jump, &
199	nx, nvex, n, ifac, la, ierr )
200	ELSE
201	CALL rpassm( work(1), work(la+1), a(ia), a(ia+ifaclainc), trigs, 1, inc, nx, &
202	jump, nvex, n, ifac, la, ierr )
203	ENDIF
204	!
205	!-- Following messages shouldn't appear in PALM applications
206	IF ( ierr /= 0 ) THEN
207	SELECT CASE (ierr)
208	CASE (:-1)
209	WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
210	CASE (0)
211	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
212	CASE (1:)
213	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
214	END SELECT
215	RETURN
216	ENDIF
217
218	la = ifac * la
219	igo = -igo
220	ia = istart
221
222	ENDDO
223	!
224	!-- If necessary, copy results back to a
225	IF ( MOD(nfax,2) /= 0 ) THEN
226	ibase = 1
227	jbase = ia
228	DO jj = 1, nvex
229	i = ibase
230	j = jbase
231	DO ii = 1, n
232	a(j) = work(i)
233	i = i + 1
234	j = j + inc
235	ENDDO
236	ibase = ibase + nx
237	jbase = jbase + jump
238	ENDDO
239	ENDIF
240	!
241	!-- Fill in zeros at end
242	ix = istart + n*inc
243	!DIR$ IVDEP
244	!CDIR NODEP
245	!OCL NOVREC
246	DO j = 1, nvex
247	a(ix) = 0.0_wp
248	a(ix+inc) = 0.0_wp
249	ix = ix + jump
250	ENDDO
251	istart = istart + nvex*jump
252	nvex = nfft
253
254	ENDDO
255
256	ELSEIF ( isign == -1 ) THEN
257	!
258	!-- Forward fft: gridpoint to spectral transform
259	istart = 1
260	DO nb = 1, nblox
261	ia = istart
262	la = n
263	igo = + 1
264
265	DO k = 1, nfax
266
267	ifac = ifax(nfax+2-k)
268	la = la / ifac
269	ierr = -1
270	IF ( igo /= -1 ) THEN
271	CALL qpassm( a(ia), a(ia+ifaclainc), work(1), work(la+1), trigs, inc, 1, jump, &
272	nx, nvex, n, ifac, la, ierr )
273	ELSE
274	CALL qpassm( work(1), work(ifacla+1), a(ia), a(ia+lainc), trigs, 1, inc, nx, &
275	jump, nvex, n, ifac, la, ierr )
276	ENDIF
277	!
278	!-- Following messages shouldn't appear in PALM applications
279	IF ( ierr /= 0 ) THEN
280	SELECT CASE (ierr)
281	CASE (:-1)
282	WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
283	CASE (0)
284	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
285	CASE (1:)
286	WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
287	END SELECT
288	RETURN
289	ENDIF
290
291	igo = -igo
292	ia = istart + inc
293
294	ENDDO
295	!
296	!-- If necessary, copy results back to a
297	IF ( MOD(nfax,2) /= 0 ) THEN
298	ibase = 1
299	jbase = ia
300	DO jj = 1, nvex
301	i = ibase
302	j = jbase
303	DO ii = 1, n
304	a(j) = work(i)
305	i = i + 1
306	j = j + inc
307	ENDDO
308	ibase = ibase + nx
309	jbase = jbase + jump
310	ENDDO
311	ENDIF
312	!
313	!-- Shift a(0) and fill in zero imag parts
314	ix = istart
315	!DIR$ IVDEP
316	!CDIR NODEP
317	!OCL NOVREC
318	DO j = 1, nvex
319	a(ix) = a(ix+inc)
320	a(ix+inc) = 0.0_wp
321	ix = ix + jump
322	ENDDO
323	IF ( MOD(n,2) /= 1 ) THEN
324	iz = istart + (n+1) * inc
325	DO j = 1, nvex
326	a(iz) = 0.0_wp
327	iz = iz + jump
328	ENDDO
329	ENDIF
330	istart = istart + nvex * jump
331	nvex = nfft
332
333	ENDDO
334
335	ENDIF
336
337	END SUBROUTINE fft991cy
338
339	!------------------------------------------------------------------------------!
340	! Description:
341	! ------------
342	!> Performs one pass through data as part of
343	!> multiple real fft (fourier analysis) routine.
344	!>
345	!> Method:
346	!>
347	!> a is first real input vector
348	!> equivalence b(1) with a(ifaclainc1+1)
349	!> c is first real output vector
350	!> equivalence d(1) with c(la*inc2+1)
351	!> trigs is a precalculated list of sines & cosines
352	!> inc1 is the addressing increment for a
353	!> inc2 is the addressing increment for c
354	!> inc3 is the increment between input vectors a
355	!> inc4 is the increment between output vectors c
356	!> lot is the number of vectors
357	!> n is the length of the vectors
358	!> ifac is the current factor of n
359	!> la = n/(product of factors used so far)
360	!> ierr is an error indicator:
361	!> 0 - pass completed without error
362	!> 1 - lot greater than nfft
363	!> 2 - ifac not catered for
364	!> 3 - ifac only catered for if la=n/ifac
365	!------------------------------------------------------------------------------!
366	SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
367
368	USE kinds
369
370	IMPLICIT NONE
371
372	INTEGER(iwp) :: ierr !<
373	INTEGER(iwp) :: ifac !<
374	INTEGER(iwp) :: inc1 !<
375	INTEGER(iwp) :: inc2 !<
376	INTEGER(iwp) :: inc3 !<
377	INTEGER(iwp) :: inc4 !<
378	INTEGER(iwp) :: la !<
379	INTEGER(iwp) :: lot !<
380	INTEGER(iwp) :: n !<
381
382	!
383	!-- Arrays are dimensioned with n
384	REAL(wp) :: a(*) !<
385	REAL(wp) :: b(*) !<
386	REAL(wp) :: c(*) !<
387	REAL(wp) :: d(*) !<
388	REAL(wp) :: trigs(*) !<
389
390	REAL(wp) :: a0 !<
391	REAL(wp) :: a1 !<
392	REAL(wp) :: a10 !<
393	REAL(wp) :: a11 !<
394	REAL(wp) :: a2 !<
395	REAL(wp) :: a20 !<
396	REAL(wp) :: a21 !<
397	REAL(wp) :: a3 !<
398	REAL(wp) :: a4 !<
399	REAL(wp) :: a5 !<
400	REAL(wp) :: a6 !<
401	REAL(wp) :: b0 !<
402	REAL(wp) :: b1 !<
403	REAL(wp) :: b10 !<
404	REAL(wp) :: b11 !<
405	REAL(wp) :: b2 !<
406	REAL(wp) :: b20 !<
407	REAL(wp) :: b21 !<
408	REAL(wp) :: b3 !<
409	REAL(wp) :: b4 !<
410	REAL(wp) :: b5 !<
411	REAL(wp) :: b6 !<
412	REAL(wp) :: c1 !<
413	REAL(wp) :: c2 !<
414	REAL(wp) :: c3 !<
415	REAL(wp) :: c4 !<
416	REAL(wp) :: c5 !<
417	REAL(wp) :: qrt5 !<
418	REAL(wp) :: s1 !<
419	REAL(wp) :: s2 !<
420	REAL(wp) :: s3 !<
421	REAL(wp) :: s4 !<
422	REAL(wp) :: s5 !<
423	REAL(wp) :: sin36 !<
424	REAL(wp) :: sin45 !<
425	REAL(wp) :: sin60 !<
426	REAL(wp) :: sin72 !<
427	REAL(wp) :: z !<
428	REAL(wp) :: zqrt5 !<
429	REAL(wp) :: zsin36 !<
430	REAL(wp) :: zsin45 !<
431	REAL(wp) :: zsin60 !<
432	REAL(wp) :: zsin72 !<
433
434	INTEGER(iwp) :: i !<
435	INTEGER(iwp) :: ia !<
436	INTEGER(iwp) :: ib !<
437	INTEGER(iwp) :: ibase !<
438	INTEGER(iwp) :: ic !<
439	INTEGER(iwp) :: id !<
440	INTEGER(iwp) :: ie !<
441	INTEGER(iwp) :: if !<
442	INTEGER(iwp) :: ig !<
443	INTEGER(iwp) :: igo !<
444	INTEGER(iwp) :: ih !<
445	INTEGER(iwp) :: iink !<
446	INTEGER(iwp) :: ijk !<
447	INTEGER(iwp) :: ijump !<
448	INTEGER(iwp) :: j !<
449	INTEGER(iwp) :: ja !<
450	INTEGER(iwp) :: jb !<
451	INTEGER(iwp) :: jbase !<
452	INTEGER(iwp) :: jc !<
453	INTEGER(iwp) :: jd !<
454	INTEGER(iwp) :: je !<
455	INTEGER(iwp) :: jf !<
456	INTEGER(iwp) :: jink !<
457	INTEGER(iwp) :: k !<
458	INTEGER(iwp) :: kb !<
459	INTEGER(iwp) :: kc !<
460	INTEGER(iwp) :: kd !<
461	INTEGER(iwp) :: ke !<
462	INTEGER(iwp) :: kf !<
463	INTEGER(iwp) :: kstop !<
464	INTEGER(iwp) :: l !<
465	INTEGER(iwp) :: m !<
466
467
468	DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
469	qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
470
471
472	ierr = 0
473
474	IF ( lot > nfft ) THEN
475	ierr = 1
476	RETURN
477	ENDIF
478
479	m = n / ifac
480	iink = la * inc1
481	jink = la * inc2
482	ijump = (ifac-1) * iink
483	kstop = ( n-ifac ) / ( 2*ifac )
484
485	ibase = 0
486	jbase = 0
487	igo = ifac - 1
488	IF ( igo == 7 ) igo = 6
489
490	IF (igo < 1 .OR. igo > 6 ) THEN
491	ierr = 2
492	RETURN
493	ENDIF
494
495
496	SELECT CASE ( igo )
497
498	!
499	!-- Coding for factor 2
500	CASE ( 1 )
501	ia = 1
502	ib = ia + iink
503	ja = 1
504	jb = ja + (2m-la) inc2
505
506	IF ( la /= m ) THEN
507
508	DO l = 1, la
509	i = ibase
510	j = jbase
511	!DIR$ IVDEP
512	!CDIR NODEP
513	!OCL NOVREC
514	DO ijk = 1, lot
515	c(ja+j) = a(ia+i) + a(ib+i)
516	c(jb+j) = a(ia+i) - a(ib+i)
517	i = i + inc3
518	j = j + inc4
519	ENDDO
520	ibase = ibase + inc1
521	jbase = jbase + inc2
522	ENDDO
523
524	ja = ja + jink
525	jink = 2 * jink
526	jb = jb - jink
527	ibase = ibase + ijump
528	ijump = 2 * ijump + iink
529
530	IF ( ja /= jb ) THEN
531
532	DO k = la, kstop, la
533	kb = k + k
534	c1 = trigs(kb+1)
535	s1 = trigs(kb+2)
536	jbase = 0
537	DO l = 1, la
538	i = ibase
539	j = jbase
540	!DIR$ IVDEP
541	!CDIR NODEP
542	!OCL NOVREC
543	DO ijk = 1, lot
544	c(ja+j) = a(ia+i) + (c1a(ib+i)+s1b(ib+i))
545	c(jb+j) = a(ia+i) - (c1a(ib+i)+s1b(ib+i))
546	d(ja+j) = (c1b(ib+i)-s1a(ib+i)) + b(ia+i)
547	d(jb+j) = (c1b(ib+i)-s1a(ib+i)) - b(ia+i)
548	i = i + inc3
549	j = j + inc4
550	ENDDO
551	ibase = ibase + inc1
552	jbase = jbase + inc2
553	ENDDO
554
555	ibase = ibase + ijump
556	ja = ja + jink
557	jb = jb - jink
558	ENDDO
559
560	IF ( ja > jb ) RETURN
561
562	ENDIF
563
564	jbase = 0
565	DO l = 1, la
566	i = ibase
567	j = jbase
568	!DIR$ IVDEP
569	!CDIR NODEP
570	!OCL NOVREC
571	DO ijk = 1, lot
572	c(ja+j) = a(ia+i)
573	d(ja+j) = -a(ib+i)
574	i = i + inc3
575	j = j + inc4
576
577	ENDDO
578	ibase = ibase + inc1
579	jbase = jbase + inc2
580	ENDDO
581
582	ELSE
583
584	z = 1.0_wp/REAL(n,KIND=wp)
585	DO l = 1, la
586	i = ibase
587	j = jbase
588	!DIR$ IVDEP
589	!CDIR NODEP
590	!OCL NOVREC
591	DO ijk = 1, lot
592	c(ja+j) = z*(a(ia+i)+a(ib+i))
593	c(jb+j) = z*(a(ia+i)-a(ib+i))
594	i = i + inc3
595	j = j + inc4
596	ENDDO
597	ibase = ibase + inc1
598	jbase = jbase + inc2
599	ENDDO
600
601	ENDIF
602
603	!
604	!-- Coding for factor 3
605	CASE ( 2 )
606
607	ia = 1
608	ib = ia + iink
609	ic = ib + iink
610	ja = 1
611	jb = ja + (2m-la) inc2
612	jc = jb
613
614	IF ( la /= m ) THEN
615
616	DO l = 1, la
617	i = ibase
618	j = jbase
619	!DIR$ IVDEP
620	!CDIR NODEP
621	!OCL NOVREC
622	DO ijk = 1, lot
623	c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
624	c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
625	d(jb+j) = sin60*(a(ic+i)-a(ib+i))
626	i = i + inc3
627	j = j + inc4
628	ENDDO
629	ibase = ibase + inc1
630	jbase = jbase + inc2
631	ENDDO
632
633	ja = ja + jink
634	jink = 2 * jink
635	jb = jb + jink
636	jc = jc - jink
637	ibase = ibase + ijump
638	ijump = 2 * ijump + iink
639
640	IF ( ja /= jc ) THEN
641
642	DO k = la, kstop, la
643	kb = k + k
644	kc = kb + kb
645	c1 = trigs(kb+1)
646	s1 = trigs(kb+2)
647	c2 = trigs(kc+1)
648	s2 = trigs(kc+2)
649
650	jbase = 0
651	DO l = 1, la
652	i = ibase
653	j = jbase
654	!DIR$ IVDEP
655	!CDIR NODEP
656	!OCL NOVREC
657	DO ijk = 1, lot
658	a1 = (c1a(ib+i)+s1b(ib+i)) + (c2a(ic+i)+s2b(ic+i))
659	b1 = (c1b(ib+i)-s1a(ib+i)) + (c2b(ic+i)-s2a(ic+i))
660	a2 = a(ia+i) - 0.5_wp*a1
661	b2 = b(ia+i) - 0.5_wp*b1
662	a3 = sin60((c1a(ib+i)+s1b(ib+i))-(c2a(ic+i)+s2*b(ic+i)))
663	b3 = sin60((c1b(ib+i)-s1a(ib+i))-(c2b(ic+i)-s2*a(ic+i)))
664	c(ja+j) = a(ia+i) + a1
665	d(ja+j) = b(ia+i) + b1
666	c(jb+j) = a2 + b3
667	d(jb+j) = b2 - a3
668	c(jc+j) = a2 - b3
669	d(jc+j) = -(b2+a3)
670	i = i + inc3
671	j = j + inc4
672	ENDDO
673	ibase = ibase + inc1
674	jbase = jbase + inc2
675	ENDDO
676	ibase = ibase + ijump
677	ja = ja + jink
678	jb = jb + jink
679	jc = jc - jink
680	ENDDO
681
682	IF ( ja > jc ) RETURN
683
684	ENDIF
685
686	jbase = 0
687	DO l = 1, la
688	i = ibase
689	j = jbase
690	!DIR$ IVDEP
691	!CDIR NODEP
692	!OCL NOVREC
693	DO ijk = 1, lot
694	c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
695	d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
696	c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
697	i = i + inc3
698	j = j + inc4
699	ENDDO
700	ibase = ibase + inc1
701	jbase = jbase + inc2
702	ENDDO
703
704	ELSE
705
706	z = 1.0_wp / REAL( n, KIND=wp )
707	zsin60 = z*sin60
708	DO l = 1, la
709	i = ibase
710	j = jbase
711	!DIR$ IVDEP
712	!CDIR NODEP
713	!OCL NOVREC
714	DO ijk = 1, lot
715	c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
716	c(jb+j) = z(a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))
717	d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
718	i = i + inc3
719	j = j + inc4
720	ENDDO
721	ibase = ibase + inc1
722	jbase = jbase + inc2
723	ENDDO
724
725	ENDIF
726
727	!
728	!-- Coding for factor 4
729	CASE ( 3 )
730
731	ia = 1
732	ib = ia + iink
733	ic = ib + iink
734	id = ic + iink
735	ja = 1
736	jb = ja + (2m-la) inc2
737	jc = jb + 2minc2
738	jd = jb
739
740	IF ( la /= m ) THEN
741
742	DO l = 1, la
743
744	i = ibase
745	j = jbase
746	!DIR$ IVDEP
747	!CDIR NODEP
748	!OCL NOVREC
749	DO ijk = 1, lot
750	c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
751	c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
752	c(jb+j) = a(ia+i) - a(ic+i)
753	d(jb+j) = a(id+i) - a(ib+i)
754	i = i + inc3
755	j = j + inc4
756	ENDDO
757	ibase = ibase + inc1
758	jbase = jbase + inc2
759	ENDDO
760
761	ja = ja + jink
762	jink = 2 * jink
763	jb = jb + jink
764	jc = jc - jink
765	jd = jd - jink
766	ibase = ibase + ijump
767	ijump = 2 * ijump + iink
768
769	IF ( jb /= jc ) THEN
770
771	DO k = la, kstop, la
772	kb = k + k
773	kc = kb + kb
774	kd = kc + kb
775	c1 = trigs(kb+1)
776	s1 = trigs(kb+2)
777	c2 = trigs(kc+1)
778	s2 = trigs(kc+2)
779	c3 = trigs(kd+1)
780	s3 = trigs(kd+2)
781	jbase = 0
782	DO l = 1, la
783	i = ibase
784	j = jbase
785	!DIR$ IVDEP
786	!CDIR NODEP
787	!OCL NOVREC
788	DO ijk = 1, lot
789	a0 = a(ia+i) + (c2a(ic+i)+s2b(ic+i))
790	a2 = a(ia+i) - (c2a(ic+i)+s2b(ic+i))
791	a1 = (c1a(ib+i)+s1b(ib+i)) + (c3a(id+i)+s3b(id+i))
792	a3 = (c1a(ib+i)+s1b(ib+i)) - (c3a(id+i)+s3b(id+i))
793	b0 = b(ia+i) + (c2b(ic+i)-s2a(ic+i))
794	b2 = b(ia+i) - (c2b(ic+i)-s2a(ic+i))
795	b1 = (c1b(ib+i)-s1a(ib+i)) + (c3b(id+i)-s3a(id+i))
796	b3 = (c1b(ib+i)-s1a(ib+i)) - (c3b(id+i)-s3a(id+i))
797	c(ja+j) = a0 + a1
798	c(jc+j) = a0 - a1
799	d(ja+j) = b0 + b1
800	d(jc+j) = b1 - b0
801	c(jb+j) = a2 + b3
802	c(jd+j) = a2 - b3
803	d(jb+j) = b2 - a3
804	d(jd+j) = -(b2+a3)
805	i = i + inc3
806	j = j + inc4
807	ENDDO
808	ibase = ibase + inc1
809	jbase = jbase + inc2
810	ENDDO
811	ibase = ibase + ijump
812	ja = ja + jink
813	jb = jb + jink
814	jc = jc - jink
815	jd = jd - jink
816	ENDDO
817
818	IF ( jb > jc ) RETURN
819
820	ENDIF
821
822	sin45 = SQRT( 0.5_wp )
823	jbase = 0
824	DO l = 1, la
825	i = ibase
826	j = jbase
827	!DIR$ IVDEP
828	!CDIR NODEP
829	!OCL NOVREC
830	DO ijk = 1, lot
831	c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
832	c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
833	d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
834	d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
835	i = i + inc3
836	j = j + inc4
837	ENDDO
838	ibase = ibase + inc1
839	jbase = jbase + inc2
840	ENDDO
841
842	ELSE
843
844	z = 1.0_wp / REAL( n, KIND=wp )
845	DO l = 1, la
846	i = ibase
847	j = jbase
848	!DIR$ IVDEP
849	!CDIR NODEP
850	!OCL NOVREC
851	DO ijk = 1, lot
852	c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
853	c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
854	c(jb+j) = z*(a(ia+i)-a(ic+i))
855	d(jb+j) = z*(a(id+i)-a(ib+i))
856	i = i + inc3
857	j = j + inc4
858	ENDDO
859	ibase = ibase + inc1
860	jbase = jbase + inc2
861	ENDDO
862
863	ENDIF
864
865	!
866	!-- Coding for factor 5
867	CASE ( 4 )
868
869	ia = 1
870	ib = ia + iink
871	ic = ib + iink
872	id = ic + iink
873	ie = id + iink
874	ja = 1
875	jb = ja + (2m-la) inc2
876	jc = jb + 2minc2
877	jd = jc
878	je = jb
879
880	IF ( la /= m ) THEN
881
882	DO l = 1, la
883	i = ibase
884	j = jbase
885	!DIR$ IVDEP
886	!CDIR NODEP
887	!OCL NOVREC
888	DO ijk = 1, lot
889	a1 = a(ib+i) + a(ie+i)
890	a3 = a(ib+i) - a(ie+i)
891	a2 = a(ic+i) + a(id+i)
892	a4 = a(ic+i) - a(id+i)
893	a5 = a(ia+i) - 0.25_wp*(a1+a2)
894	a6 = qrt5*(a1-a2)
895	c(ja+j) = a(ia+i) + (a1+a2)
896	c(jb+j) = a5 + a6
897	c(jc+j) = a5 - a6
898	d(jb+j) = -sin72a3 - sin36a4
899	d(jc+j) = -sin36a3 + sin72a4
900	i = i + inc3
901	j = j + inc4
902	ENDDO
903	ibase = ibase + inc1
904	jbase = jbase + inc2
905	ENDDO
906
907	ja = ja + jink
908	jink = 2 * jink
909	jb = jb + jink
910	jc = jc + jink
911	jd = jd - jink
912	je = je - jink
913	ibase = ibase + ijump
914	ijump = 2 * ijump + iink
915
916	IF ( jb /= jd ) THEN
917
918	DO k = la, kstop, la
919	kb = k + k
920	kc = kb + kb
921	kd = kc + kb
922	ke = kd + kb
923	c1 = trigs(kb+1)
924	s1 = trigs(kb+2)
925	c2 = trigs(kc+1)
926	s2 = trigs(kc+2)
927	c3 = trigs(kd+1)
928	s3 = trigs(kd+2)
929	c4 = trigs(ke+1)
930	s4 = trigs(ke+2)
931	jbase = 0
932	DO l = 1, la
933	i = ibase
934	j = jbase
935	!DIR$ IVDEP
936	!CDIR NODEP
937	!OCL NOVREC
938	DO ijk = 1, lot
939	a1 = (c1a(ib+i)+s1b(ib+i)) + (c4a(ie+i)+s4b(ie+i))
940	a3 = (c1a(ib+i)+s1b(ib+i)) - (c4a(ie+i)+s4b(ie+i))
941	a2 = (c2a(ic+i)+s2b(ic+i)) + (c3a(id+i)+s3b(id+i))
942	a4 = (c2a(ic+i)+s2b(ic+i)) - (c3a(id+i)+s3b(id+i))
943	b1 = (c1b(ib+i)-s1a(ib+i)) + (c4b(ie+i)-s4a(ie+i))
944	b3 = (c1b(ib+i)-s1a(ib+i)) - (c4b(ie+i)-s4a(ie+i))
945	b2 = (c2b(ic+i)-s2a(ic+i)) + (c3b(id+i)-s3a(id+i))
946	b4 = (c2b(ic+i)-s2a(ic+i)) - (c3b(id+i)-s3a(id+i))
947	a5 = a(ia+i) - 0.25_wp*(a1+a2)
948	a6 = qrt5*(a1-a2)
949	b5 = b(ia+i) - 0.25_wp*(b1+b2)
950	b6 = qrt5*(b1-b2)
951	a10 = a5 + a6
952	a20 = a5 - a6
953	b10 = b5 + b6
954	b20 = b5 - b6
955	a11 = sin72b3 + sin36b4
956	a21 = sin36b3 - sin72b4
957	b11 = sin72a3 + sin36a4
958	b21 = sin36a3 - sin72a4
959	c(ja+j) = a(ia+i) + (a1+a2)
960	c(jb+j) = a10 + a11
961	c(je+j) = a10 - a11
962	c(jc+j) = a20 + a21
963	c(jd+j) = a20 - a21
964	d(ja+j) = b(ia+i) + (b1+b2)
965	d(jb+j) = b10 - b11
966	d(je+j) = -(b10+b11)
967	d(jc+j) = b20 - b21
968	d(jd+j) = -(b20+b21)
969	i = i + inc3
970	j = j + inc4
971	ENDDO
972	ibase = ibase + inc1
973	jbase = jbase + inc2
974	ENDDO
975	ibase = ibase + ijump
976	ja = ja + jink
977	jb = jb + jink
978	jc = jc + jink
979	jd = jd - jink
980	je = je - jink
981	ENDDO
982
983	IF ( jb > jd ) RETURN
984
985	ENDIF
986
987	jbase = 0
988	DO l = 1, la
989	i = ibase
990	j = jbase
991	!DIR$ IVDEP
992	!CDIR NODEP
993	!OCL NOVREC
994	DO ijk = 1, lot
995	a1 = a(ib+i) + a(ie+i)
996	a3 = a(ib+i) - a(ie+i)
997	a2 = a(ic+i) + a(id+i)
998	a4 = a(ic+i) - a(id+i)
999	a5 = a(ia+i) + 0.25_wp*(a3-a4)
1000	a6 = qrt5*(a3+a4)
1001	c(ja+j) = a5 + a6
1002	c(jb+j) = a5 - a6
1003	c(jc+j) = a(ia+i) - (a3-a4)
1004	d(ja+j) = -sin36a1 - sin72a2
1005	d(jb+j) = -sin72a1 + sin36a2
1006	i = i + inc3
1007	j = j + inc4
1008	ENDDO
1009	ibase = ibase + inc1
1010	jbase = jbase + inc2
1011	ENDDO
1012
1013	ELSE
1014
1015	z = 1.0_wp / REAL( n, KIND=wp )
1016	zqrt5 = z * qrt5
1017	zsin36 = z * sin36
1018	zsin72 = z * sin72
1019	DO l = 1, la
1020	i = ibase
1021	j = jbase
1022	!DIR$ IVDEP
1023	!CDIR NODEP
1024	!OCL NOVREC
1025	DO ijk = 1, lot
1026	a1 = a(ib+i) + a(ie+i)
1027	a3 = a(ib+i) - a(ie+i)
1028	a2 = a(ic+i) + a(id+i)
1029	a4 = a(ic+i) - a(id+i)
1030	a5 = z(a(ia+i)-0.25_wp(a1+a2))
1031	a6 = zqrt5*(a1-a2)
1032	c(ja+j) = z*(a(ia+i)+(a1+a2))
1033	c(jb+j) = a5 + a6
1034	c(jc+j) = a5 - a6
1035	d(jb+j) = -zsin72a3 - zsin36a4
1036	d(jc+j) = -zsin36a3 + zsin72a4
1037	i = i + inc3
1038	j = j + inc4
1039	ENDDO
1040	ibase = ibase + inc1
1041	jbase = jbase + inc2
1042	ENDDO
1043
1044	ENDIF
1045
1046	!
1047	!-- Coding for factor 6
1048	CASE ( 5 )
1049
1050	ia = 1
1051	ib = ia + iink
1052	ic = ib + iink
1053	id = ic + iink
1054	ie = id + iink
1055	if = ie + iink
1056	ja = 1
1057	jb = ja + (2m-la) inc2
1058	jc = jb + 2minc2
1059	jd = jc + 2minc2
1060	je = jc
1061	jf = jb
1062
1063	IF ( la /= m ) THEN
1064
1065	DO l = 1, la
1066	i = ibase
1067	j = jbase
1068	!DIR$ IVDEP
1069	!CDIR NODEP
1070	!OCL NOVREC
1071	DO ijk = 1, lot
1072	a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
1073	c(ja+j) = (a(ia+i)+a(id+i)) + a11
1074	c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
1075	d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
1076	a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
1077	c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
1078	d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1079	c(jd+j) = (a(ia+i)-a(id+i)) + a11
1080	i = i + inc3
1081	j = j + inc4
1082	ENDDO
1083	ibase = ibase + inc1
1084	jbase = jbase + inc2
1085	ENDDO
1086	ja = ja + jink
1087	jink = 2 * jink
1088	jb = jb + jink
1089	jc = jc + jink
1090	jd = jd - jink
1091	je = je - jink
1092	jf = jf - jink
1093	ibase = ibase + ijump
1094	ijump = 2 * ijump + iink
1095
1096	IF ( jc /= jd ) THEN
1097
1098	DO k = la, kstop, la
1099	kb = k + k
1100	kc = kb + kb
1101	kd = kc + kb
1102	ke = kd + kb
1103	kf = ke + kb
1104	c1 = trigs(kb+1)
1105	s1 = trigs(kb+2)
1106	c2 = trigs(kc+1)
1107	s2 = trigs(kc+2)
1108	c3 = trigs(kd+1)
1109	s3 = trigs(kd+2)
1110	c4 = trigs(ke+1)
1111	s4 = trigs(ke+2)
1112	c5 = trigs(kf+1)
1113	s5 = trigs(kf+2)
1114	jbase = 0
1115	DO l = 1, la
1116	i = ibase
1117	j = jbase
1118	!DIR$ IVDEP
1119	!CDIR NODEP
1120	!OCL NOVREC
1121	DO ijk = 1, lot
1122	a1 = c1a(ib+i) + s1b(ib+i)
1123	b1 = c1b(ib+i) - s1a(ib+i)
1124	a2 = c2a(ic+i) + s2b(ic+i)
1125	b2 = c2b(ic+i) - s2a(ic+i)
1126	a3 = c3a(id+i) + s3b(id+i)
1127	b3 = c3b(id+i) - s3a(id+i)
1128	a4 = c4a(ie+i) + s4b(ie+i)
1129	b4 = c4b(ie+i) - s4a(ie+i)
1130	a5 = c5a(if+i) + s5b(if+i)
1131	b5 = c5b(if+i) - s5a(if+i)
1132	a11 = (a2+a5) + (a1+a4)
1133	a20 = (a(ia+i)+a3) - 0.5_wp*a11
1134	a21 = sin60*((a2+a5)-(a1+a4))
1135	b11 = (b2+b5) + (b1+b4)
1136	b20 = (b(ia+i)+b3) - 0.5_wp*b11
1137	b21 = sin60*((b2+b5)-(b1+b4))
1138	c(ja+j) = (a(ia+i)+a3) + a11
1139	d(ja+j) = (b(ia+i)+b3) + b11
1140	c(jc+j) = a20 - b21
1141	d(jc+j) = a21 + b20
1142	c(je+j) = a20 + b21
1143	d(je+j) = a21 - b20
1144	a11 = (a2-a5) + (a4-a1)
1145	a20 = (a(ia+i)-a3) - 0.5_wp*a11
1146	a21 = sin60*((a4-a1)-(a2-a5))
1147	b11 = (b5-b2) - (b4-b1)
1148	b20 = (b3-b(ia+i)) - 0.5_wp*b11
1149	b21 = sin60*((b5-b2)+(b4-b1))
1150	c(jb+j) = a20 - b21
1151	d(jb+j) = a21 - b20
1152	c(jd+j) = a11 + (a(ia+i)-a3)
1153	d(jd+j) = b11 + (b3-b(ia+i))
1154	c(jf+j) = a20 + b21
1155	d(jf+j) = a21 + b20
1156	i = i + inc3
1157	j = j + inc4
1158	ENDDO
1159	ibase = ibase + inc1
1160	jbase = jbase + inc2
1161	ENDDO
1162	ibase = ibase + ijump
1163	ja = ja + jink
1164	jb = jb + jink
1165	jc = jc + jink
1166	jd = jd - jink
1167	je = je - jink
1168	jf = jf - jink
1169	ENDDO
1170
1171	IF ( jc > jd ) RETURN
1172
1173	ENDIF
1174
1175	jbase = 0
1176	DO l = 1, la
1177	i = ibase
1178	j = jbase
1179	!DIR$ IVDEP
1180	!CDIR NODEP
1181	!OCL NOVREC
1182	DO ijk = 1, lot
1183	c(ja+j) = (a(ia+i)+0.5_wp(a(ic+i)-a(ie+i))) + sin60(a(ib+i)-a(if+i))
1184	d(ja+j) = -(a(id+i)+0.5_wp(a(ib+i)+a(if+i))) - sin60(a(ic+i)+a(ie+i))
1185	c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
1186	d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
1187	c(jc+j) = (a(ia+i)+0.5_wp(a(ic+i)-a(ie+i))) - sin60(a(ib+i)-a(if+i))
1188	d(jc+j) = -(a(id+i)+0.5_wp(a(ib+i)+a(if+i))) + sin60(a(ic+i)+a(ie+i))
1189	i = i + inc3
1190	j = j + inc4
1191	ENDDO
1192	ibase = ibase + inc1
1193	jbase = jbase + inc2
1194	ENDDO
1195
1196	ELSE
1197
1198	z = 1.0_wp/REAL(n,KIND=wp)
1199	zsin60 = z*sin60
1200	DO l = 1, la
1201	i = ibase
1202	j = jbase
1203	!DIR$ IVDEP
1204	!CDIR NODEP
1205	!OCL NOVREC
1206	DO ijk = 1, lot
1207	a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
1208	c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
1209	c(jc+j) = z((a(ia+i)+a(id+i))-0.5_wpa11)
1210	d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
1211	a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
1212	c(jb+j) = z((a(ia+i)-a(id+i))-0.5_wpa11)
1213	d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
1214	c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
1215	i = i + inc3
1216	j = j + inc4
1217	ENDDO
1218	ibase = ibase + inc1
1219	jbase = jbase + inc2
1220	ENDDO
1221
1222	ENDIF
1223
1224	!
1225	!-- Coding for factor 8
1226	CASE ( 6 )
1227
1228	IF ( la /= m ) THEN
1229	ierr = 3
1230	RETURN
1231	ENDIF
1232
1233	ia = 1
1234	ib = ia + iink
1235	ic = ib + iink
1236	id = ic + iink
1237	ie = id + iink
1238	if = ie + iink
1239	ig = if + iink
1240	ih = ig + iink
1241	ja = 1
1242	jb = ja + la * inc2
1243	jc = jb + 2minc2
1244	jd = jc + 2minc2
1245	je = jd + 2minc2
1246	z = 1.0_wp / REAL( n, KIND=wp )
1247	zsin45 = z * SQRT( 0.5_wp )
1248
1249	DO l = 1, la
1250	i = ibase
1251	j = jbase
1252
1253	!DIR$ IVDEP
1254	!CDIR NODEP
1255	!OCL NOVREC
1256	DO ijk = 1, lot
1257	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))))
1258	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))))
1259	c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
1260	d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
1261	c(jb+j) = z(a(ia+i)-a(ie+i)) + zsin45((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
1262	c(jd+j) = z(a(ia+i)-a(ie+i)) - zsin45((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
1263	d(jb+j) = zsin45((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z(a(ig+i)-a(ic+i))
1264	d(jd+j) = zsin45((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z(a(ig+i)-a(ic+i))
1265	i = i + inc3
1266	j = j + inc4
1267	ENDDO
1268	ibase = ibase + inc1
1269	jbase = jbase + inc2
1270	ENDDO
1271
1272	END SELECT
1273
1274	END SUBROUTINE qpassm
1275
1276	!------------------------------------------------------------------------------!
1277	! Description:
1278	! ------------
1279	!> Same as qpassm, but for backward fft
1280	!------------------------------------------------------------------------------!
1281	SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
1282
1283	USE kinds
1284
1285	IMPLICIT NONE
1286
1287	INTEGER(iwp) :: ierr !<
1288	INTEGER(iwp) :: ifac !<
1289	INTEGER(iwp) :: inc1 !<
1290	INTEGER(iwp) :: inc2 !<
1291	INTEGER(iwp) :: inc3 !<
1292	INTEGER(iwp) :: inc4 !<
1293	INTEGER(iwp) :: la !<
1294	INTEGER(iwp) :: lot !<
1295	INTEGER(iwp) :: n !<
1296	!
1297	!-- Arrays are dimensioned with n
1298	REAL(wp) :: a(*) !<
1299	REAL(wp) :: b(*) !<
1300	REAL(wp) :: c(*) !<
1301	REAL(wp) :: d(*) !<
1302	REAL(wp) :: trigs(*) !<
1303
1304	REAL(wp) :: c1 !<
1305	REAL(wp) :: c2 !<
1306	REAL(wp) :: c3 !<
1307	REAL(wp) :: c4 !<
1308	REAL(wp) :: c5 !<
1309	REAL(wp) :: qqrt5 !<
1310	REAL(wp) :: qrt5 !<
1311	REAL(wp) :: s1 !<
1312	REAL(wp) :: s2 !<
1313	REAL(wp) :: s3 !<
1314	REAL(wp) :: s4 !<
1315	REAL(wp) :: s5 !<
1316	REAL(wp) :: sin36 !<
1317	REAL(wp) :: sin45 !<
1318	REAL(wp) :: sin60 !<
1319	REAL(wp) :: sin72 !<
1320	REAL(wp) :: ssin36 !<
1321	REAL(wp) :: ssin45 !<
1322	REAL(wp) :: ssin60 !<
1323	REAL(wp) :: ssin72 !<
1324
1325	INTEGER(iwp) :: i !<
1326	INTEGER(iwp) :: ia !<
1327	INTEGER(iwp) :: ib !<
1328	INTEGER(iwp) :: ibase !<
1329	INTEGER(iwp) :: ic !<
1330	INTEGER(iwp) :: id !<
1331	INTEGER(iwp) :: ie !<
1332	INTEGER(iwp) :: if !<
1333	INTEGER(iwp) :: igo !<
1334	INTEGER(iwp) :: iink !<
1335	INTEGER(iwp) :: ijk !<
1336	INTEGER(iwp) :: j !<
1337	INTEGER(iwp) :: ja !<
1338	INTEGER(iwp) :: jb !<
1339	INTEGER(iwp) :: jbase !<
1340	INTEGER(iwp) :: jc !<
1341	INTEGER(iwp) :: jd !<
1342	INTEGER(iwp) :: je !<
1343	INTEGER(iwp) :: jf !<
1344	INTEGER(iwp) :: jg !<
1345	INTEGER(iwp) :: jh !<
1346	INTEGER(iwp) :: jink !<
1347	INTEGER(iwp) :: jump !<
1348	INTEGER(iwp) :: k !<
1349	INTEGER(iwp) :: kb !<
1350	INTEGER(iwp) :: kc !<
1351	INTEGER(iwp) :: kd !<
1352	INTEGER(iwp) :: ke !<
1353	INTEGER(iwp) :: kf !<
1354	INTEGER(iwp) :: kstop !<
1355	INTEGER(iwp) :: l !<
1356	INTEGER(iwp) :: m !<
1357
1358	REAL(wp) :: a10(nfft) !<
1359	REAL(wp) :: a11(nfft) !<
1360	REAL(wp) :: a20(nfft) !<
1361	REAL(wp) :: a21(nfft) !<
1362	REAL(wp) :: b10(nfft) !<
1363	REAL(wp) :: b11(nfft) !<
1364	REAL(wp) :: b20(nfft) !<
1365	REAL(wp) :: b21(nfft) !<
1366
1367
1368	DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
1369	qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
1370
1371
1372	ierr = 0
1373
1374	m = n / ifac
1375	iink = la * inc1
1376	jink = la * inc2
1377	jump = (ifac-1) * jink
1378	kstop = (n-ifac) / (2*ifac)
1379
1380	IF ( lot > nfft ) THEN
1381	ierr = 1
1382	RETURN
1383	ENDIF
1384	ibase = 0
1385	jbase = 0
1386	igo = ifac - 1
1387	IF ( igo == 7 ) igo = 6
1388	IF ( igo < 1 .OR. igo > 6 ) THEN
1389	ierr = 2
1390	RETURN
1391	ENDIF
1392
1393
1394	SELECT CASE ( igo )
1395	!
1396	!-- Coding for factor 2
1397	CASE ( 1 )
1398
1399	ia = 1
1400	ib = ia + (2m-la) inc1
1401	ja = 1
1402	jb = ja + jink
1403
1404	IF ( la /= m ) THEN
1405
1406	DO l = 1, la
1407	i = ibase
1408	j = jbase
1409	!DIR$ IVDEP
1410	!CDIR NODEP
1411	!OCL NOVREC
1412	DO ijk = 1, lot
1413	c(ja+j) = a(ia+i) + a(ib+i)
1414	c(jb+j) = a(ia+i) - a(ib+i)
1415	i = i + inc3
1416	j = j + inc4
1417	ENDDO
1418	ibase = ibase + inc1
1419	jbase = jbase + inc2
1420	ENDDO
1421
1422	ia = ia + iink
1423	iink = 2 * iink
1424	ib = ib - iink
1425	ibase = 0
1426	jbase = jbase + jump
1427	jump = 2 * jump + jink
1428
1429	IF ( ia /= ib ) THEN
1430
1431	DO k = la, kstop, la
1432	kb = k + k
1433	c1 = trigs(kb+1)
1434	s1 = trigs(kb+2)
1435	ibase = 0
1436	DO l = 1, la
1437	i = ibase
1438	j = jbase
1439	!DIR$ IVDEP
1440	!CDIR NODEP
1441	!OCL NOVREC
1442	DO ijk = 1, lot
1443	c(ja+j) = a(ia+i) + a(ib+i)
1444	d(ja+j) = b(ia+i) - b(ib+i)
1445	c(jb+j) = c1(a(ia+i)-a(ib+i)) - s1(b(ia+i)+b(ib+i))
1446	d(jb+j) = s1(a(ia+i)-a(ib+i)) + c1(b(ia+i)+b(ib+i))
1447	i = i + inc3
1448	j = j + inc4
1449	ENDDO
1450	ibase = ibase + inc1
1451	jbase = jbase + inc2
1452	ENDDO
1453	ia = ia + iink
1454	ib = ib - iink
1455	jbase = jbase + jump
1456	ENDDO
1457
1458	IF ( ia > ib ) RETURN
1459
1460	ENDIF
1461
1462	ibase = 0
1463	DO l = 1, la
1464	i = ibase
1465	j = jbase
1466	!DIR$ IVDEP
1467	!CDIR NODEP
1468	!OCL NOVREC
1469	DO ijk = 1, lot
1470	c(ja+j) = a(ia+i)
1471	c(jb+j) = -b(ia+i)
1472	i = i + inc3
1473	j = j + inc4
1474	ENDDO
1475	ibase = ibase + inc1
1476	jbase = jbase + inc2
1477	ENDDO
1478
1479	ELSE
1480
1481	DO l = 1, la
1482	i = ibase
1483	j = jbase
1484	!DIR$ IVDEP
1485	!CDIR NODEP
1486	!OCL NOVREC
1487	DO ijk = 1, lot
1488	c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
1489	c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
1490	i = i + inc3
1491	j = j + inc4
1492	ENDDO
1493	ibase = ibase + inc1
1494	jbase = jbase + inc2
1495	ENDDO
1496
1497	ENDIF
1498
1499	!
1500	!-- Coding for factor 3
1501	CASE ( 2 )
1502
1503	ia = 1
1504	ib = ia + (2m-la) inc1
1505	ic = ib
1506	ja = 1
1507	jb = ja + jink
1508	jc = jb + jink
1509
1510	IF ( la /= m ) THEN
1511
1512	DO l = 1, la
1513	i = ibase
1514	j = jbase
1515	!DIR$ IVDEP
1516	!CDIR NODEP
1517	!OCL NOVREC
1518	DO ijk = 1, lot
1519	c(ja+j) = a(ia+i) + a(ib+i)
1520	c(jb+j) = (a(ia+i)-0.5_wpa(ib+i)) - (sin60(b(ib+i)))
1521	c(jc+j) = (a(ia+i)-0.5_wpa(ib+i)) + (sin60(b(ib+i)))
1522	i = i + inc3
1523	j = j + inc4
1524	ENDDO
1525	ibase = ibase + inc1
1526	jbase = jbase + inc2
1527	ENDDO
1528	ia = ia + iink
1529	iink = 2 * iink
1530	ib = ib + iink
1531	ic = ic - iink
1532	jbase = jbase + jump
1533	jump = 2 * jump + jink
1534
1535	IF ( ia /= ic ) THEN
1536
1537	DO k = la, kstop, la
1538	kb = k + k
1539	kc = kb + kb
1540	c1 = trigs(kb+1)
1541	s1 = trigs(kb+2)
1542	c2 = trigs(kc+1)
1543	s2 = trigs(kc+2)
1544	ibase = 0
1545	DO l = 1, la
1546	i = ibase
1547	j = jbase
1548	!DIR$ IVDEP
1549	!CDIR NODEP
1550	!OCL NOVREC
1551	DO ijk = 1, lot
1552	c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1553	d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
1554	c(jb+j) = c1((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
1555	- s1((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
1556	d(jb+j) = s1((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
1557	+ c1((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
1558	c(jc+j) = c2((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
1559	- s2((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
1560	d(jc+j) = s2((a(ia+i)-0.5_wp(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
1561	+ c2((b(ia+i)-0.5_wp(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
1562	i = i + inc3
1563	j = j + inc4
1564	ENDDO
1565	ibase = ibase + inc1
1566	jbase = jbase + inc2
1567	ENDDO
1568	ia = ia + iink
1569	ib = ib + iink
1570	ic = ic - iink
1571	jbase = jbase + jump
1572	ENDDO
1573
1574	IF ( ia > ic ) RETURN
1575
1576	ENDIF
1577
1578	ibase = 0
1579	DO l = 1, la
1580	i = ibase
1581	j = jbase
1582	!DIR$ IVDEP
1583	!CDIR NODEP
1584	!OCL NOVREC
1585	DO ijk = 1, lot
1586	c(ja+j) = a(ia+i) + a(ib+i)
1587	c(jb+j) = (0.5_wpa(ia+i)-a(ib+i)) - (sin60b(ia+i))
1588	c(jc+j) = -(0.5_wpa(ia+i)-a(ib+i)) - (sin60b(ia+i))
1589	i = i + inc3
1590	j = j + inc4
1591	ENDDO
1592	ibase = ibase + inc1
1593	jbase = jbase + inc2
1594	ENDDO
1595
1596	ELSE
1597
1598	ssin60 = 2.0_wp * sin60
1599	DO l = 1, la
1600	i = ibase
1601	j = jbase
1602	!DIR$ IVDEP
1603	!CDIR NODEP
1604	!OCL NOVREC
1605	DO ijk = 1, lot
1606	c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
1607	c(jb+j) = (2.0_wpa(ia+i)-a(ib+i)) - (ssin60b(ib+i))
1608	c(jc+j) = (2.0_wpa(ia+i)-a(ib+i)) + (ssin60b(ib+i))
1609	i = i + inc3
1610	j = j + inc4
1611	ENDDO
1612	ibase = ibase + inc1
1613	jbase = jbase + inc2
1614	ENDDO
1615
1616	ENDIF
1617
1618	!
1619	!-- Coding for factor 4
1620	CASE ( 3 )
1621
1622	ia = 1
1623	ib = ia + (2m-la) inc1
1624	ic = ib + 2minc1
1625	id = ib
1626	ja = 1
1627	jb = ja + jink
1628	jc = jb + jink
1629	jd = jc + jink
1630
1631	IF ( la /= m ) THEN
1632
1633	DO l = 1, la
1634	i = ibase
1635	j = jbase
1636	!DIR$ IVDEP
1637	!CDIR NODEP
1638	!OCL NOVREC
1639	DO ijk = 1, lot
1640	c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
1641	c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
1642	c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
1643	c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
1644	i = i + inc3
1645	j = j + inc4
1646	ENDDO
1647	ibase = ibase + inc1
1648	jbase = jbase + inc2
1649	ENDDO
1650	ia = ia + iink
1651	iink = 2 * iink
1652	ib = ib + iink
1653	ic = ic - iink
1654	id = id - iink
1655	jbase = jbase + jump
1656	jump = 2 * jump + jink
1657
1658	IF ( ib /= ic ) THEN
1659
1660	DO k = la, kstop, la
1661	kb = k + k
1662	kc = kb + kb
1663	kd = kc + kb
1664	c1 = trigs(kb+1)
1665	s1 = trigs(kb+2)
1666	c2 = trigs(kc+1)
1667	s2 = trigs(kc+2)
1668	c3 = trigs(kd+1)
1669	s3 = trigs(kd+2)
1670	ibase = 0
1671	DO l = 1, la
1672	i = ibase
1673	j = jbase
1674	!DIR$ IVDEP
1675	!CDIR NODEP
1676	!OCL NOVREC
1677	DO ijk = 1, lot
1678	c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
1679	d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
1680	c(jc+j) = c2((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2((b(ia+i)-b(ic+i))&
1681	-(b(ib+i)-b(id+i)))
1682	d(jc+j) = s2((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2((b(ia+i)-b(ic+i))&
1683	-(b(ib+i)-b(id+i)))
1684	c(jb+j) = c1((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1((b(ia+i)+b(ic+i))&
1685	+(a(ib+i)-a(id+i)))
1686	d(jb+j) = s1((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1((b(ia+i)+b(ic+i))&
1687	+(a(ib+i)-a(id+i)))
1688	c(jd+j) = c3((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3((b(ia+i)+b(ic+i))&
1689	-(a(ib+i)-a(id+i)))
1690	d(jd+j) = s3((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3((b(ia+i)+b(ic+i))&
1691	-(a(ib+i)-a(id+i)))
1692	i = i + inc3
1693	j = j + inc4
1694	ENDDO
1695	ibase = ibase + inc1
1696	jbase = jbase + inc2
1697	ENDDO
1698	ia = ia + iink
1699	ib = ib + iink
1700	ic = ic - iink
1701	id = id - iink
1702	jbase = jbase + jump
1703	ENDDO
1704
1705	IF ( ib > ic ) RETURN
1706
1707	ENDIF
1708
1709	ibase = 0
1710	sin45 = SQRT( 0.5_wp )
1711	DO l = 1, la
1712	i = ibase
1713	j = jbase
1714	!DIR$ IVDEP
1715	!CDIR NODEP
1716	!OCL NOVREC
1717	DO ijk = 1, lot
1718	c(ja+j) = a(ia+i) + a(ib+i)
1719	c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
1720	c(jc+j) = b(ib+i) - b(ia+i)
1721	c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
1722	i = i + inc3
1723	j = j + inc4
1724	ENDDO
1725	ibase = ibase + inc1
1726	jbase = jbase + inc2
1727	ENDDO
1728
1729	ELSE
1730
1731	DO l = 1, la
1732	i = ibase
1733	j = jbase
1734	!DIR$ IVDEP
1735	!CDIR NODEP
1736	!OCL NOVREC
1737	DO ijk = 1, lot
1738	c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
1739	c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
1740	c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
1741	c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
1742	i = i + inc3
1743	j = j + inc4
1744	ENDDO
1745	ibase = ibase + inc1
1746	jbase = jbase + inc2
1747	ENDDO
1748
1749	ENDIF
1750
1751	!
1752	!-- Coding for factor 5
1753	CASE ( 4 )
1754
1755	ia = 1
1756	ib = ia + (2m-la) inc1
1757	ic = ib + 2minc1
1758	id = ic
1759	ie = ib
1760	ja = 1
1761	jb = ja + jink
1762	jc = jb + jink
1763	jd = jc + jink
1764	je = jd + jink
1765
1766	IF ( la /= m ) THEN
1767
1768	DO l = 1, la
1769	i = ibase
1770	j = jbase
1771	!DIR$ IVDEP
1772	!CDIR NODEP
1773	!OCL NOVREC
1774	DO ijk = 1, lot
1775	c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
1776	c(jb+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))+qrt5(a(ib+i)-a(ic+i))) - &
1777	(sin72b(ib+i)+sin36b(ic+i))
1778	c(jc+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))-qrt5(a(ib+i)-a(ic+i))) - &
1779	(sin36b(ib+i)-sin72b(ic+i))
1780	c(jd+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))-qrt5(a(ib+i)-a(ic+i))) + &
1781	(sin36b(ib+i)-sin72b(ic+i))
1782	c(je+j) = ((a(ia+i)-0.25_wp(a(ib+i)+a(ic+i)))+qrt5(a(ib+i)-a(ic+i))) + &
1783	(sin72b(ib+i)+sin36b(ic+i))
1784	i = i + inc3
1785	j = j + inc4
1786	ENDDO
1787	ibase = ibase + inc1
1788	jbase = jbase + inc2
1789	ENDDO
1790	ia = ia + iink
1791	iink = 2 * iink
1792	ib = ib + iink
1793	ic = ic + iink
1794	id = id - iink
1795	ie = ie - iink
1796	jbase = jbase + jump
1797	jump = 2 * jump + jink
1798
1799	IF ( ib /= id ) THEN
1800
1801	DO k = la, kstop, la
1802	kb = k + k
1803	kc = kb + kb
1804	kd = kc + kb
1805	ke = kd + kb
1806	c1 = trigs(kb+1)
1807	s1 = trigs(kb+2)
1808	c2 = trigs(kc+1)
1809	s2 = trigs(kc+2)
1810	c3 = trigs(kd+1)
1811	s3 = trigs(kd+2)
1812	c4 = trigs(ke+1)
1813	s4 = trigs(ke+2)
1814	ibase = 0
1815	DO l = 1, la
1816	i = ibase
1817	j = jbase
1818	!DIR$ IVDEP
1819	!CDIR NODEP
1820	!OCL NOVREC
1821	DO ijk = 1, lot
1822	a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
1823	qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1824	a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
1825	qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
1826	b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
1827	qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1828	b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
1829	qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
1830	a11(ijk) = sin72(b(ib+i)+b(ie+i)) + sin36(b(ic+i)+b(id+i))
1831	a21(ijk) = sin36(b(ib+i)+b(ie+i)) - sin72(b(ic+i)+b(id+i))
1832	b11(ijk) = sin72(a(ib+i)-a(ie+i)) + sin36(a(ic+i)-a(id+i))
1833	b21(ijk) = sin36(a(ib+i)-a(ie+i)) - sin72(a(ic+i)-a(id+i))
1834	c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
1835	d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
1836	c(jb+j) = c1(a10(ijk)-a11(ijk)) - s1(b10(ijk)+b11(ijk))
1837	d(jb+j) = s1(a10(ijk)-a11(ijk)) + c1(b10(ijk)+b11(ijk))
1838	c(je+j) = c4(a10(ijk)+a11(ijk)) - s4(b10(ijk)-b11(ijk))
1839	d(je+j) = s4(a10(ijk)+a11(ijk)) + c4(b10(ijk)-b11(ijk))
1840	c(jc+j) = c2(a20(ijk)-a21(ijk)) - s2(b20(ijk)+b21(ijk))
1841	d(jc+j) = s2(a20(ijk)-a21(ijk)) + c2(b20(ijk)+b21(ijk))
1842	c(jd+j) = c3(a20(ijk)+a21(ijk)) - s3(b20(ijk)-b21(ijk))
1843	d(jd+j) = s3(a20(ijk)+a21(ijk)) + c3(b20(ijk)-b21(ijk))
1844
1845	i = i + inc3
1846	j = j + inc4
1847	ENDDO
1848	ibase = ibase + inc1
1849	jbase = jbase + inc2
1850	ENDDO
1851	ia = ia + iink
1852	ib = ib + iink
1853	ic = ic + iink
1854	id = id - iink
1855	ie = ie - iink
1856	jbase = jbase + jump
1857	ENDDO
1858
1859	IF ( ib > id ) RETURN
1860
1861	ENDIF
1862
1863	ibase = 0
1864	DO l = 1, la
1865	i = ibase
1866	j = jbase
1867	!DIR$ IVDEP
1868	!CDIR NODEP
1869	!OCL NOVREC
1870	DO ijk = 1, lot
1871	c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
1872	c(jb+j) = (qrt5(a(ia+i)-a(ib+i))+(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
1873	(sin36b(ia+i)+sin72b(ib+i))
1874	c(je+j) = -(qrt5(a(ia+i)-a(ib+i))+(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
1875	(sin36b(ia+i)+sin72b(ib+i))
1876	c(jc+j) = (qrt5(a(ia+i)-a(ib+i))-(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
1877	(sin72b(ia+i)-sin36b(ib+i))
1878	c(jd+j) = -(qrt5(a(ia+i)-a(ib+i))-(0.25_wp(a(ia+i)+a(ib+i))-a(ic+i))) - &
1879	(sin72b(ia+i)-sin36b(ib+i))
1880	i = i + inc3
1881	j = j + inc4
1882	ENDDO
1883	ibase = ibase + inc1
1884	jbase = jbase + inc2
1885	ENDDO
1886
1887	ELSE
1888
1889	qqrt5 = 2.0_wp * qrt5
1890	ssin36 = 2.0_wp * sin36
1891	ssin72 = 2.0_wp * sin72
1892	DO l = 1, la
1893	i = ibase
1894	j = jbase
1895	!DIR$ IVDEP
1896	!CDIR NODEP
1897	!OCL NOVREC
1898	DO ijk = 1, lot
1899	c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
1900	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))) &
1901	- (ssin72b(ib+i)+ssin36b(ic+i))
1902	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))) &
1903	- (ssin36b(ib+i)-ssin72b(ic+i))
1904	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))) &
1905	+ (ssin36b(ib+i)-ssin72b(ic+i))
1906	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))) &
1907	+ (ssin72b(ib+i)+ssin36b(ic+i))
1908	i = i + inc3
1909	j = j + inc4
1910	ENDDO
1911	ibase = ibase + inc1
1912	jbase = jbase + inc2
1913	ENDDO
1914
1915	ENDIF
1916
1917	!
1918	!-- Coding for factor 6
1919	CASE ( 5 )
1920
1921	ia = 1
1922	ib = ia + (2m-la) inc1
1923	ic = ib + 2minc1
1924	id = ic + 2minc1
1925	ie = ic
1926	if = ib
1927	ja = 1
1928	jb = ja + jink
1929	jc = jb + jink
1930	jd = jc + jink
1931	je = jd + jink
1932	jf = je + jink
1933
1934	IF ( la /= m ) THEN
1935
1936	DO l = 1, la
1937	i = ibase
1938	j = jbase
1939	!DIR$ IVDEP
1940	!CDIR NODEP
1941	!OCL NOVREC
1942	DO ijk = 1, lot
1943	c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
1944	c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
1945	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)))
1946	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)))
1947	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)))
1948	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)))
1949	i = i + inc3
1950	j = j + inc4
1951	ENDDO
1952	ibase = ibase + inc1
1953	jbase = jbase + inc2
1954	ENDDO
1955	ia = ia + iink
1956	iink = 2 * iink
1957	ib = ib + iink
1958	ic = ic + iink
1959	id = id - iink
1960	ie = ie - iink
1961	if = if - iink
1962	jbase = jbase + jump
1963	jump = 2 * jump + jink
1964
1965	IF ( ic /= id ) THEN
1966
1967	DO k = la, kstop, la
1968	kb = k + k
1969	kc = kb + kb
1970	kd = kc + kb
1971	ke = kd + kb
1972	kf = ke + kb
1973	c1 = trigs(kb+1)
1974	s1 = trigs(kb+2)
1975	c2 = trigs(kc+1)
1976	s2 = trigs(kc+2)
1977	c3 = trigs(kd+1)
1978	s3 = trigs(kd+2)
1979	c4 = trigs(ke+1)
1980	s4 = trigs(ke+2)
1981	c5 = trigs(kf+1)
1982	s5 = trigs(kf+2)
1983	ibase = 0
1984	DO l = 1, la
1985	i = ibase
1986	j = jbase
1987	!DIR$ IVDEP
1988	!CDIR NODEP
1989	!OCL NOVREC
1990	DO ijk = 1, lot
1991	a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(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	b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
1995	b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
1996	b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
1997
1998	c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
1999	d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
2000	c(jc+j) = c2(a20(ijk)-b21(ijk)) - s2(b20(ijk)+a21(ijk))
2001	d(jc+j) = s2(a20(ijk)-b21(ijk)) + c2(b20(ijk)+a21(ijk))
2002	c(je+j) = c4(a20(ijk)+b21(ijk)) - s4(b20(ijk)-a21(ijk))
2003	d(je+j) = s4(a20(ijk)+b21(ijk)) + c4(b20(ijk)-a21(ijk))
2004
2005	a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
2006	b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
2007	a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
2008	a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
2009	b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
2010	b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
2011
2012	c(jd+j) = c3((a(ia+i)-a(id+i))+a11(ijk)) - s3((b(ia+i)+b(id+i))-b11(ijk))
2013	d(jd+j) = s3((a(ia+i)-a(id+i))+a11(ijk)) + c3((b(ia+i)+b(id+i))-b11(ijk))
2014	c(jb+j) = c1(a20(ijk)-b21(ijk)) - s1(b20(ijk)-a21(ijk))
2015	d(jb+j) = s1(a20(ijk)-b21(ijk)) + c1(b20(ijk)-a21(ijk))
2016	c(jf+j) = c5(a20(ijk)+b21(ijk)) - s5(b20(ijk)+a21(ijk))
2017	d(jf+j) = s5(a20(ijk)+b21(ijk)) + c5(b20(ijk)+a21(ijk))
2018
2019	i = i + inc3
2020	j = j + inc4
2021	ENDDO
2022	ibase = ibase + inc1
2023	jbase = jbase + inc2
2024	ENDDO
2025	ia = ia + iink
2026	ib = ib + iink
2027	ic = ic + iink
2028	id = id - iink
2029	ie = ie - iink
2030	if = if - iink
2031	jbase = jbase + jump
2032	ENDDO
2033
2034	IF ( ic > id ) RETURN
2035
2036	ENDIF
2037
2038	ibase = 0
2039	DO l = 1, la
2040	i = ibase
2041	j = jbase
2042	!DIR$ IVDEP
2043	!CDIR NODEP
2044	!OCL NOVREC
2045	DO ijk = 1, lot
2046	c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
2047	c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
2048	c(jb+j) = (sin60(a(ia+i)-a(ic+i))) - (0.5_wp(b(ia+i)+b(ic+i))+b(ib+i))
2049	c(jf+j) = -(sin60(a(ia+i)-a(ic+i))) - (0.5_wp(b(ia+i)+b(ic+i))+b(ib+i))
2050	c(jc+j) = sin60(b(ic+i)-b(ia+i)) + (0.5_wp(a(ia+i)+a(ic+i))-a(ib+i))
2051	c(je+j) = sin60(b(ic+i)-b(ia+i)) - (0.5_wp(a(ia+i)+a(ic+i))-a(ib+i))
2052	i = i + inc3
2053	j = j + inc4
2054	ENDDO
2055	ibase = ibase + inc1
2056	jbase = jbase + inc2
2057	ENDDO
2058
2059	ELSE
2060
2061	ssin60 = 2.0_wp * sin60
2062	DO l = 1, la
2063	i = ibase
2064	j = jbase
2065	!DIR$ IVDEP
2066	!CDIR NODEP
2067	!OCL NOVREC
2068	DO ijk = 1, lot
2069	c(ja+j) = (2.0_wp(a(ia+i)+a(id+i))) + (2.0_wp(a(ib+i)+a(ic+i)))
2070	c(jd+j) = (2.0_wp(a(ia+i)-a(id+i))) - (2.0_wp(a(ib+i)-a(ic+i)))
2071	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)))
2072	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)))
2073	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)))
2074	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)))
2075	i = i + inc3
2076	j = j + inc4
2077	ENDDO
2078	ibase = ibase + inc1
2079	jbase = jbase + inc2
2080	ENDDO
2081
2082	ENDIF
2083
2084	!
2085	!-- Coding for factor 8
2086	CASE ( 6 )
2087
2088	IF ( la /= m ) THEN
2089	ierr = 3
2090	RETURN
2091	ENDIF
2092
2093	ia = 1
2094	ib = ia + la*inc1
2095	ic = ib + 2lainc1
2096	id = ic + 2lainc1
2097	ie = id + 2lainc1
2098	ja = 1
2099	jb = ja + jink
2100	jc = jb + jink
2101	jd = jc + jink
2102	je = jd + jink
2103	jf = je + jink
2104	jg = jf + jink
2105	jh = jg + jink
2106	ssin45 = SQRT( 2.0_wp )
2107
2108	DO l = 1, la
2109	i = ibase
2110	j = jbase
2111	!DIR$ IVDEP
2112	!CDIR NODEP
2113	!OCL NOVREC
2114	DO ijk = 1, lot
2115	c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
2116	c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
2117	c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
2118	c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
2119	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)))
2120	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)))
2121	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)))
2122	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)))
2123	i = i + inc3
2124	j = j + inc4
2125	ENDDO
2126	ibase = ibase + inc1
2127	jbase = jbase + inc2
2128	ENDDO
2129
2130	END SELECT
2131
2132	END SUBROUTINE rpassm
2133
2134	!------------------------------------------------------------------------------!
2135	! Description:
2136	! ------------
2137	!> Computes factors of n & trigonometric functins required by fft99 & fft991cy
2138	!> Method: Dimension trigs(n),ifax(1),jfax(10),lfax(7)
2139	!> subroutine 'set99' - computes factors of n & trigonometric
2140	!> functins required by fft99 & fft991cy
2141	!------------------------------------------------------------------------------!
2142	SUBROUTINE set99( trigs, ifax, n )
2143
2144
2145	USE control_parameters, &
2146	ONLY: message_string
2147
2148	USE kinds
2149
2150	IMPLICIT NONE
2151
2152	INTEGER(iwp) :: n !<
2153
2154	INTEGER(iwp) :: ifax(*) !<
2155	REAL(wp) :: trigs(*) !<
2156
2157
2158	REAL(wp) :: angle !<
2159	REAL(wp) :: del !<
2160	INTEGER(iwp) :: i !<
2161	INTEGER(iwp) :: ifac !<
2162	INTEGER(iwp) :: ixxx !<
2163	INTEGER(iwp) :: k !<
2164	INTEGER(iwp) :: l !<
2165	INTEGER(iwp) :: nfax !<
2166	INTEGER(iwp) :: nhl !<
2167	INTEGER(iwp) :: nil !<
2168	INTEGER(iwp) :: nu !<
2169
2170	INTEGER(iwp) :: jfax(10) !<
2171	INTEGER(iwp) :: lfax(7) !<
2172
2173	DATA lfax / 6, 8, 5, 4, 3, 2, 1 /
2174
2175
2176	ixxx = 1
2177
2178	del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp )
2179	nil = 0
2180	nhl = (n/2) - 1
2181	DO k = nil, nhl
2182	angle = REAL( k, KIND=wp ) * del
2183	trigs(2*k+1) = COS( angle )
2184	trigs(2*k+2) = SIN( angle )
2185	ENDDO
2186
2187	!
2188	!-- Find factors of n (8,6,5,4,3,2; only one 8 allowed).
2189	!-- Look for sixes first, store factors in descending order.
2190	nu = n
2191	ifac = 6
2192	k = 0
2193	l = 1
2194
2195	DO
2196
2197	IF ( MOD( nu, ifac ) == 0 ) THEN
2198	k = k + 1
2199	jfax(k) = ifac
2200	IF ( ifac == 8 .AND. k /= 1 ) THEN
2201	jfax(1) = 8
2202	jfax(k) = 6
2203	ENDIF
2204	nu = nu/ifac
2205	IF ( nu == 1 ) EXIT
2206	IF ( ifac /= 8 ) CYCLE
2207	ENDIF
2208
2209	l = l + 1
2210	ifac = lfax(l)
2211	IF (ifac < 2 ) THEN
2212	message_string = 'number of gridpoints along x or/and y ' // &
2213	'contain illegal factors' // &
2214	'&only factors 2, 3, 5 are allowed'
2215	CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
2216	RETURN
2217	ENDIF
2218
2219	ENDDO
2220	!
2221	!-- Now reverse order of factors
2222	nfax = k
2223	ifax(1) = nfax
2224	DO i = 1, nfax
2225	ifax(nfax+2-i) = jfax(i)
2226	ENDDO
2227	ifax(10) = n
2228
2229	END SUBROUTINE set99
2230
2231	END MODULE temperton_fft

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |