Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy.f90 @ 1191

Last change on this file since 1191 was 1167, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 46.0 KB

Rev	Line
[1]	1	MODULE fft_xy
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
	22	!
[1167]	23	!
[1]	24	! Former revisions:
	25	! -----------------
[3]	26	! $Id: fft_xy.f90 1167 2013-05-24 14:03:34Z heinze $
[392]	27	!
[1167]	28	! 1166 2013-05-24 13:55:44Z raasch
	29	! C_DOUBLE/COMPLEX reset to dpk
	30	!
[1154]	31	! 1153 2013-05-10 14:33:08Z raasch
	32	! code adjustment of data types for CUDA fft required by PGI 12.3 / CUDA 5.0
	33	!
[1112]	34	! 1111 2013-03-08 23:54:10Z raasch
	35	! further openACC statements added, CUDA branch completely runs on GPU
	36	! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
	37	! used total domain)
	38	!
[1107]	39	! 1106 2013-03-04 05:31:38Z raasch
	40	! CUDA fft added
	41	! array_kind renamed precision_kind, 3D- instead of 1D-loops in fft_x and fft_y
	42	! old fft_x, fft_y become fft_x_1d, fft_y_1d and are used for 1D-decomposition
	43	!
[1093]	44	! 1092 2013-02-02 11:24:22Z raasch
	45	! variable sizw declared for NEC case only
	46	!
[1037]	47	! 1036 2012-10-22 13:43:42Z raasch
	48	! code put under GPL (PALM 3.9)
	49	!
[392]	50	! 274 2009-03-26 15:11:21Z heinze
	51	! Output of messages replaced by message handling routine.
	52	!
	53	! Feb. 2007
[3]	54	! RCS Log replace by Id keyword, revision history cleaned up
	55	!
[1]	56	! Revision 1.4 2006/03/28 12:27:09 raasch
	57	! Stop when system-specific fft is selected on NEC. For unknown reasons this
	58	! causes a program abort during first allocation in init_grid.
	59	!
	60	! Revision 1.2 2004/04/30 11:44:27 raasch
	61	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
	62	! fft_x and fft_y,
	63	! function FFT replaced by subroutine FFTN due to problems with 64-bit
	64	! mode on ibm,
	65	! shape of array cwork is explicitly stored in ishape/jshape and handled
	66	! to routine FFTN instead of shape-function (due to compiler error on
	67	! decalpha),
	68	! non vectorized FFT for nec included
	69	!
	70	! Revision 1.1 2002/06/11 13:00:49 raasch
	71	! Initial revision
	72	!
	73	!
	74	! Description:
	75	! ------------
	76	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
	77	! Original version: Klaus Ketelsen (May 2002)
	78	!------------------------------------------------------------------------------!
	79
	80	USE control_parameters
	81	USE indices
[1153]	82	#if defined( __cuda_fft )
	83	USE ISO_C_BINDING
	84	#endif
[1106]	85	USE precision_kind
[1]	86	USE singleton
	87	USE temperton_fft
[1106]	88	USE transpose_indices
[1]	89
	90	IMPLICIT NONE
	91
	92	PRIVATE
[1106]	93	PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
[1]	94
	95	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
	96
	97	LOGICAL, SAVE :: init_fft = .FALSE.
	98
[1106]	99	REAL, SAVE :: dnx, dny, sqr_dnx, sqr_dny
[1]	100	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
	101
	102	#if defined( __ibm )
	103	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
	104	!
	105	!-- The following working arrays contain tables and have to be "save" and
	106	!-- shared in OpenMP sense
	107	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
	108	#elif defined( __nec )
	109	INTEGER, SAVE :: nz1
	110	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
	111	trig_yf
[1106]	112	#elif defined( __cuda_fft )
[1153]	113	INTEGER(C_INT), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
	114	INTEGER, SAVE :: total_points_x_transpo, total_points_y_transpo
[1]	115	#endif
	116
	117	!
	118	!-- Public interfaces
	119	INTERFACE fft_init
	120	MODULE PROCEDURE fft_init
	121	END INTERFACE fft_init
	122
	123	INTERFACE fft_x
	124	MODULE PROCEDURE fft_x
	125	END INTERFACE fft_x
	126
[1106]	127	INTERFACE fft_x_1d
	128	MODULE PROCEDURE fft_x_1d
	129	END INTERFACE fft_x_1d
	130
[1]	131	INTERFACE fft_y
	132	MODULE PROCEDURE fft_y
	133	END INTERFACE fft_y
	134
[1106]	135	INTERFACE fft_y_1d
	136	MODULE PROCEDURE fft_y_1d
	137	END INTERFACE fft_y_1d
	138
[1]	139	INTERFACE fft_x_m
	140	MODULE PROCEDURE fft_x_m
	141	END INTERFACE fft_x_m
	142
	143	INTERFACE fft_y_m
	144	MODULE PROCEDURE fft_y_m
	145	END INTERFACE fft_y_m
	146
	147	CONTAINS
	148
	149
	150	SUBROUTINE fft_init
	151
[1106]	152	USE cuda_fft_interfaces
	153
[1]	154	IMPLICIT NONE
	155
	156	!
	157	!-- The following temporary working arrays have to be on stack or private
	158	!-- in OpenMP sense
	159	#if defined( __ibm )
	160	REAL, DIMENSION(0:nx+2) :: workx
	161	REAL, DIMENSION(0:ny+2) :: worky
	162	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
	163	#elif defined( __nec )
	164	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
	165	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
	166	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
	167	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
	168	#endif
	169
	170	!
	171	!-- Return, if already called
	172	IF ( init_fft ) THEN
	173	RETURN
	174	ELSE
	175	init_fft = .TRUE.
	176	ENDIF
	177
	178	IF ( fft_method == 'system-specific' ) THEN
	179
[1106]	180	dnx = 1.0 / ( nx + 1.0 )
	181	dny = 1.0 / ( ny + 1.0 )
	182	sqr_dnx = SQRT( dnx )
	183	sqr_dny = SQRT( dny )
[1]	184	#if defined( __ibm ) && ! defined( __ibmy_special )
	185	!
	186	!-- Initialize tables for fft along x
[1106]	187	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
[1]	188	aux2, nau2 )
[1106]	189	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]	190	aux4, nau2 )
	191	!
	192	!-- Initialize tables for fft along y
[1106]	193	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
[1]	194	auy2, nau2 )
[1106]	195	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]	196	auy4, nau2 )
	197	#elif defined( __nec )
[254]	198	message_string = 'fft method "' // TRIM( fft_method) // &
	199	'" currently does not work on NEC'
	200	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
[1]	201
	202	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
	203	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
	204
	205	work_x = 0.0
	206	work_y = 0.0
	207	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
	208	! when using the NEC ffts
	209
	210	!
	211	!-- Initialize tables for fft along x (non-vector and vector case (M))
[1106]	212	CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
	213	CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
	214	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
[1]	215	trig_xf, workx, 0 )
[1106]	216	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
[1]	217	trig_xb, workx, 0 )
	218	!
	219	!-- Initialize tables for fft along y (non-vector and vector case (M))
[1106]	220	CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
	221	CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
	222	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
[1]	223	trig_yf, worky, 0 )
[1106]	224	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
[1]	225	trig_yb, worky, 0 )
[1106]	226	#elif defined( __cuda_fft )
	227	total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
	228	total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
[1111]	229	CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
	230	CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
	231	CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
	232	CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
[1]	233	#else
[254]	234	message_string = 'no system-specific fft-call available'
	235	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
[1]	236	#endif
	237	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	238	!
	239	!-- Temperton-algorithm
	240	!-- Initialize tables for fft along x and y
	241	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
	242
	243	CALL set99( trigs_x, ifax_x, nx+1 )
	244	CALL set99( trigs_y, ifax_y, ny+1 )
	245
	246	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
	247
	248	CONTINUE
	249
	250	ELSE
	251
[254]	252	message_string = 'fft method "' // TRIM( fft_method) // &
	253	'" not available'
	254	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
[1]	255	ENDIF
	256
	257	END SUBROUTINE fft_init
	258
	259
	260	SUBROUTINE fft_x( ar, direction )
	261
	262	!----------------------------------------------------------------------!
	263	! fft_x !
	264	! !
	265	! Fourier-transformation along x-direction !
[1106]	266	! Version for 2D-decomposition !
[1]	267	! !
	268	! fft_x uses internal algorithms (Singleton or Temperton) or !
	269	! system-specific routines, if they are available !
	270	!----------------------------------------------------------------------!
	271
[1106]	272	USE cuda_fft_interfaces
[1153]	273	#if defined( __cuda_fft )
	274	USE ISO_C_BINDING
	275	#endif
[1106]	276
[1]	277	IMPLICIT NONE
	278
	279	CHARACTER (LEN=*) :: direction
[1111]	280	INTEGER :: i, ishape(1), j, k
[1106]	281
	282	LOGICAL :: forward_fft
	283
	284	REAL, DIMENSION(0:nx+2) :: work
	285	REAL, DIMENSION(nx+2) :: work1
	286	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	287	#if defined( __ibm )
	288	REAL, DIMENSION(nau2) :: aux2, aux4
	289	#elif defined( __nec )
	290	REAL, DIMENSION(6*(nx+1)) :: work2
	291	#elif defined( __cuda_fft )
[1111]	292	!$acc declare create( ar_tmp )
[1166]	293	COMPLEX(dpk), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: ar_tmp
[1106]	294	#endif
	295	REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar
	296
	297	IF ( direction == 'forward' ) THEN
	298	forward_fft = .TRUE.
	299	ELSE
	300	forward_fft = .FALSE.
	301	ENDIF
	302
	303	IF ( fft_method == 'singleton-algorithm' ) THEN
	304
	305	!
	306	!-- Performing the fft with singleton's software works on every system,
	307	!-- since it is part of the model
	308	ALLOCATE( cwork(0:nx) )
	309
	310	IF ( forward_fft ) then
	311
	312	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
	313	!$OMP DO
	314	DO k = nzb_x, nzt_x
	315	DO j = nys_x, nyn_x
	316
	317	DO i = 0, nx
	318	cwork(i) = CMPLX( ar(i,j,k) )
	319	ENDDO
	320
	321	ishape = SHAPE( cwork )
	322	CALL FFTN( cwork, ishape )
	323
	324	DO i = 0, (nx+1)/2
	325	ar(i,j,k) = REAL( cwork(i) )
	326	ENDDO
	327	DO i = 1, (nx+1)/2 - 1
	328	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
	329	ENDDO
	330
	331	ENDDO
	332	ENDDO
	333	!$OMP END PARALLEL
	334
	335	ELSE
	336
	337	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
	338	!$OMP DO
	339	DO k = nzb_x, nzt_x
	340	DO j = nys_x, nyn_x
	341
	342	cwork(0) = CMPLX( ar(0,j,k), 0.0 )
	343	DO i = 1, (nx+1)/2 - 1
	344	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
	345	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
	346	ENDDO
	347	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
	348
	349	ishape = SHAPE( cwork )
	350	CALL FFTN( cwork, ishape, inv = .TRUE. )
	351
	352	DO i = 0, nx
	353	ar(i,j,k) = REAL( cwork(i) )
	354	ENDDO
	355
	356	ENDDO
	357	ENDDO
	358	!$OMP END PARALLEL
	359
	360	ENDIF
	361
	362	DEALLOCATE( cwork )
	363
	364	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	365
	366	!
	367	!-- Performing the fft with Temperton's software works on every system,
	368	!-- since it is part of the model
	369	IF ( forward_fft ) THEN
	370
	371	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	372	!$OMP DO
	373	DO k = nzb_x, nzt_x
	374	DO j = nys_x, nyn_x
	375
	376	work(0:nx) = ar(0:nx,j,k)
	377	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	378
	379	DO i = 0, (nx+1)/2
	380	ar(i,j,k) = work(2*i)
	381	ENDDO
	382	DO i = 1, (nx+1)/2 - 1
	383	ar(nx+1-i,j,k) = work(2*i+1)
	384	ENDDO
	385
	386	ENDDO
	387	ENDDO
	388	!$OMP END PARALLEL
	389
	390	ELSE
	391
	392	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	393	!$OMP DO
	394	DO k = nzb_x, nzt_x
	395	DO j = nys_x, nyn_x
	396
	397	DO i = 0, (nx+1)/2
	398	work(2*i) = ar(i,j,k)
	399	ENDDO
	400	DO i = 1, (nx+1)/2 - 1
	401	work(2*i+1) = ar(nx+1-i,j,k)
	402	ENDDO
	403	work(1) = 0.0
	404	work(nx+2) = 0.0
	405
	406	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	407	ar(0:nx,j,k) = work(0:nx)
	408
	409	ENDDO
	410	ENDDO
	411	!$OMP END PARALLEL
	412
	413	ENDIF
	414
	415	ELSEIF ( fft_method == 'system-specific' ) THEN
	416
	417	#if defined( __ibm ) && ! defined( __ibmy_special )
	418	IF ( forward_fft ) THEN
	419
	420	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	421	!$OMP DO
	422	DO k = nzb_x, nzt_x
	423	DO j = nys_x, nyn_x
	424
	425	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
	426	aux2, nau2 )
	427
	428	DO i = 0, (nx+1)/2
	429	ar(i,j,k) = work(2*i)
	430	ENDDO
	431	DO i = 1, (nx+1)/2 - 1
	432	ar(nx+1-i,j,k) = work(2*i+1)
	433	ENDDO
	434
	435	ENDDO
	436	ENDDO
	437	!$OMP END PARALLEL
	438
	439	ELSE
	440
	441	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	442	!$OMP DO
	443	DO k = nzb_x, nzt_x
	444	DO j = nys_x, nyn_x
	445
	446	DO i = 0, (nx+1)/2
	447	work(2*i) = ar(i,j,k)
	448	ENDDO
	449	DO i = 1, (nx+1)/2 - 1
	450	work(2*i+1) = ar(nx+1-i,j,k)
	451	ENDDO
	452	work(1) = 0.0
	453	work(nx+2) = 0.0
	454
	455	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
	456	aux4, nau2 )
	457
	458	DO i = 0, nx
	459	ar(i,j,k) = work(i)
	460	ENDDO
	461
	462	ENDDO
	463	ENDDO
	464	!$OMP END PARALLEL
	465
	466	ENDIF
	467
	468	#elif defined( __nec )
	469
	470	IF ( forward_fft ) THEN
	471
	472	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	473	!$OMP DO
	474	DO k = nzb_x, nzt_x
	475	DO j = nys_x, nyn_x
	476
	477	work(0:nx) = ar(0:nx,j,k)
	478
	479	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
	480
	481	DO i = 0, (nx+1)/2
	482	ar(i,j,k) = work(2*i)
	483	ENDDO
	484	DO i = 1, (nx+1)/2 - 1
	485	ar(nx+1-i,j,k) = work(2*i+1)
	486	ENDDO
	487
	488	ENDDO
	489	ENDDO
	490	!$END OMP PARALLEL
	491
	492	ELSE
	493
	494	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	495	!$OMP DO
	496	DO k = nzb_x, nzt_x
	497	DO j = nys_x, nyn_x
	498
	499	DO i = 0, (nx+1)/2
	500	work(2*i) = ar(i,j,k)
	501	ENDDO
	502	DO i = 1, (nx+1)/2 - 1
	503	work(2*i+1) = ar(nx+1-i,j,k)
	504	ENDDO
	505	work(1) = 0.0
	506	work(nx+2) = 0.0
	507
	508	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
	509
	510	ar(0:nx,j,k) = work(0:nx)
	511
	512	ENDDO
	513	ENDDO
	514	!$OMP END PARALLEL
	515
	516	ENDIF
	517
	518	#elif defined( __cuda_fft )
	519
	520	IF ( forward_fft ) THEN
	521
[1111]	522	!$acc data present( ar )
	523	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
[1106]	524
[1111]	525	!$acc kernels
	526	!$acc loop
[1106]	527	DO k = nzb_x, nzt_x
	528	DO j = nys_x, nyn_x
	529
[1111]	530	!$acc loop vector( 32 )
[1106]	531	DO i = 0, (nx+1)/2
[1111]	532	ar(i,j,k) = REAL( ar_tmp(i,j,k) ) * dnx
[1106]	533	ENDDO
	534
[1111]	535	!$acc loop vector( 32 )
[1106]	536	DO i = 1, (nx+1)/2 - 1
[1111]	537	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
[1106]	538	ENDDO
	539
	540	ENDDO
	541	ENDDO
[1111]	542	!$acc end kernels
	543	!$acc end data
[1106]	544
	545	ELSE
	546
[1111]	547	!$acc data present( ar )
	548	!$acc kernels
	549	!$acc loop
[1106]	550	DO k = nzb_x, nzt_x
	551	DO j = nys_x, nyn_x
	552
[1111]	553	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
[1106]	554
[1111]	555	!$acc loop vector( 32 )
[1106]	556	DO i = 1, (nx+1)/2 - 1
[1111]	557	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
[1106]	558	ENDDO
[1111]	559	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
[1106]	560
	561	ENDDO
	562	ENDDO
[1111]	563	!$acc end kernels
[1106]	564
[1111]	565	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
	566	!$acc end data
[1106]	567
	568	ENDIF
	569
	570	#else
	571	message_string = 'no system-specific fft-call available'
	572	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
	573	#endif
	574
	575	ELSE
	576
	577	message_string = 'fft method "' // TRIM( fft_method) // &
	578	'" not available'
	579	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
	580
	581	ENDIF
	582
	583	END SUBROUTINE fft_x
	584
	585	SUBROUTINE fft_x_1d( ar, direction )
	586
	587	!----------------------------------------------------------------------!
	588	! fft_x_1d !
	589	! !
	590	! Fourier-transformation along x-direction !
	591	! Version for 1D-decomposition !
	592	! !
	593	! fft_x uses internal algorithms (Singleton or Temperton) or !
	594	! system-specific routines, if they are available !
	595	!----------------------------------------------------------------------!
	596
	597	IMPLICIT NONE
	598
	599	CHARACTER (LEN=*) :: direction
[1]	600	INTEGER :: i, ishape(1)
	601
[1106]	602	LOGICAL :: forward_fft
	603
[1]	604	REAL, DIMENSION(0:nx) :: ar
	605	REAL, DIMENSION(0:nx+2) :: work
	606	REAL, DIMENSION(nx+2) :: work1
	607	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	608	#if defined( __ibm )
	609	REAL, DIMENSION(nau2) :: aux2, aux4
	610	#elif defined( __nec )
	611	REAL, DIMENSION(6*(nx+1)) :: work2
	612	#endif
	613
[1106]	614	IF ( direction == 'forward' ) THEN
	615	forward_fft = .TRUE.
	616	ELSE
	617	forward_fft = .FALSE.
	618	ENDIF
	619
[1]	620	IF ( fft_method == 'singleton-algorithm' ) THEN
	621
	622	!
	623	!-- Performing the fft with singleton's software works on every system,
	624	!-- since it is part of the model
	625	ALLOCATE( cwork(0:nx) )
	626
[1106]	627	IF ( forward_fft ) then
[1]	628
	629	DO i = 0, nx
	630	cwork(i) = CMPLX( ar(i) )
	631	ENDDO
	632	ishape = SHAPE( cwork )
	633	CALL FFTN( cwork, ishape )
	634	DO i = 0, (nx+1)/2
	635	ar(i) = REAL( cwork(i) )
	636	ENDDO
	637	DO i = 1, (nx+1)/2 - 1
	638	ar(nx+1-i) = -AIMAG( cwork(i) )
	639	ENDDO
	640
	641	ELSE
	642
	643	cwork(0) = CMPLX( ar(0), 0.0 )
	644	DO i = 1, (nx+1)/2 - 1
	645	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
	646	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
	647	ENDDO
	648	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
	649
	650	ishape = SHAPE( cwork )
	651	CALL FFTN( cwork, ishape, inv = .TRUE. )
	652
	653	DO i = 0, nx
	654	ar(i) = REAL( cwork(i) )
	655	ENDDO
	656
	657	ENDIF
	658
	659	DEALLOCATE( cwork )
	660
	661	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	662
	663	!
	664	!-- Performing the fft with Temperton's software works on every system,
	665	!-- since it is part of the model
[1106]	666	IF ( forward_fft ) THEN
[1]	667
	668	work(0:nx) = ar
	669	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	670
	671	DO i = 0, (nx+1)/2
	672	ar(i) = work(2*i)
	673	ENDDO
	674	DO i = 1, (nx+1)/2 - 1
	675	ar(nx+1-i) = work(2*i+1)
	676	ENDDO
	677
	678	ELSE
	679
	680	DO i = 0, (nx+1)/2
	681	work(2*i) = ar(i)
	682	ENDDO
	683	DO i = 1, (nx+1)/2 - 1
	684	work(2*i+1) = ar(nx+1-i)
	685	ENDDO
	686	work(1) = 0.0
	687	work(nx+2) = 0.0
	688
	689	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	690	ar = work(0:nx)
	691
	692	ENDIF
	693
	694	ELSEIF ( fft_method == 'system-specific' ) THEN
	695
	696	#if defined( __ibm ) && ! defined( __ibmy_special )
[1106]	697	IF ( forward_fft ) THEN
[1]	698
[1106]	699	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
[1]	700	aux2, nau2 )
	701
	702	DO i = 0, (nx+1)/2
	703	ar(i) = work(2*i)
	704	ENDDO
	705	DO i = 1, (nx+1)/2 - 1
	706	ar(nx+1-i) = work(2*i+1)
	707	ENDDO
	708
	709	ELSE
	710
	711	DO i = 0, (nx+1)/2
	712	work(2*i) = ar(i)
	713	ENDDO
	714	DO i = 1, (nx+1)/2 - 1
	715	work(2*i+1) = ar(nx+1-i)
	716	ENDDO
	717	work(1) = 0.0
	718	work(nx+2) = 0.0
	719
[1106]	720	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]	721	aux4, nau2 )
	722
	723	DO i = 0, nx
	724	ar(i) = work(i)
	725	ENDDO
	726
	727	ENDIF
	728	#elif defined( __nec )
[1106]	729	IF ( forward_fft ) THEN
[1]	730
	731	work(0:nx) = ar(0:nx)
	732
[1106]	733	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
	734
[1]	735	DO i = 0, (nx+1)/2
	736	ar(i) = work(2*i)
	737	ENDDO
	738	DO i = 1, (nx+1)/2 - 1
	739	ar(nx+1-i) = work(2*i+1)
	740	ENDDO
	741
	742	ELSE
	743
	744	DO i = 0, (nx+1)/2
	745	work(2*i) = ar(i)
	746	ENDDO
	747	DO i = 1, (nx+1)/2 - 1
	748	work(2*i+1) = ar(nx+1-i)
	749	ENDDO
	750	work(1) = 0.0
	751	work(nx+2) = 0.0
	752
[1106]	753	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
[1]	754
	755	ar(0:nx) = work(0:nx)
	756
	757	ENDIF
	758	#else
[254]	759	message_string = 'no system-specific fft-call available'
[1106]	760	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
[1]	761	#endif
	762	ELSE
[274]	763	message_string = 'fft method "' // TRIM( fft_method) // &
	764	'" not available'
[1106]	765	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]	766
	767	ENDIF
	768
[1106]	769	END SUBROUTINE fft_x_1d
[1]	770
	771	SUBROUTINE fft_y( ar, direction )
	772
	773	!----------------------------------------------------------------------!
	774	! fft_y !
	775	! !
	776	! Fourier-transformation along y-direction !
[1106]	777	! Version for 2D-decomposition !
[1]	778	! !
	779	! fft_y uses internal algorithms (Singleton or Temperton) or !
	780	! system-specific routines, if they are available !
	781	!----------------------------------------------------------------------!
	782
[1106]	783	USE cuda_fft_interfaces
[1153]	784	#if defined( __cuda_fft )
	785	USE ISO_C_BINDING
	786	#endif
[1106]	787
[1]	788	IMPLICIT NONE
	789
	790	CHARACTER (LEN=*) :: direction
[1111]	791	INTEGER :: i, j, jshape(1), k
[1106]	792
	793	LOGICAL :: forward_fft
	794
	795	REAL, DIMENSION(0:ny+2) :: work
	796	REAL, DIMENSION(ny+2) :: work1
	797	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	798	#if defined( __ibm )
	799	REAL, DIMENSION(nau2) :: auy2, auy4
	800	#elif defined( __nec )
	801	REAL, DIMENSION(6*(ny+1)) :: work2
	802	#elif defined( __cuda_fft )
[1111]	803	!$acc declare create( ar_tmp )
[1166]	804	COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp
[1106]	805	#endif
	806	REAL, DIMENSION(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) :: ar
	807
	808	IF ( direction == 'forward' ) THEN
	809	forward_fft = .TRUE.
	810	ELSE
	811	forward_fft = .FALSE.
	812	ENDIF
	813
	814	IF ( fft_method == 'singleton-algorithm' ) THEN
	815
	816	!
	817	!-- Performing the fft with singleton's software works on every system,
	818	!-- since it is part of the model
	819	ALLOCATE( cwork(0:ny) )
	820
	821	IF ( forward_fft ) then
	822
	823	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
	824	!$OMP DO
	825	DO k = nzb_y, nzt_y
	826	DO i = nxl_y, nxr_y
	827
	828	DO j = 0, ny
	829	cwork(j) = CMPLX( ar(j,i,k) )
	830	ENDDO
	831
	832	jshape = SHAPE( cwork )
	833	CALL FFTN( cwork, jshape )
	834
	835	DO j = 0, (ny+1)/2
	836	ar(j,i,k) = REAL( cwork(j) )
	837	ENDDO
	838	DO j = 1, (ny+1)/2 - 1
	839	ar(ny+1-j,i,k) = -AIMAG( cwork(j) )
	840	ENDDO
	841
	842	ENDDO
	843	ENDDO
	844	!$OMP END PARALLEL
	845
	846	ELSE
	847
	848	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
	849	!$OMP DO
	850	DO k = nzb_y, nzt_y
	851	DO i = nxl_y, nxr_y
	852
	853	cwork(0) = CMPLX( ar(0,i,k), 0.0 )
	854	DO j = 1, (ny+1)/2 - 1
	855	cwork(j) = CMPLX( ar(j,i,k), -ar(ny+1-j,i,k) )
	856	cwork(ny+1-j) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
	857	ENDDO
	858	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
	859
	860	jshape = SHAPE( cwork )
	861	CALL FFTN( cwork, jshape, inv = .TRUE. )
	862
	863	DO j = 0, ny
	864	ar(j,i,k) = REAL( cwork(j) )
	865	ENDDO
	866
	867	ENDDO
	868	ENDDO
	869	!$OMP END PARALLEL
	870
	871	ENDIF
	872
	873	DEALLOCATE( cwork )
	874
	875	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	876
	877	!
	878	!-- Performing the fft with Temperton's software works on every system,
	879	!-- since it is part of the model
	880	IF ( forward_fft ) THEN
	881
	882	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	883	!$OMP DO
	884	DO k = nzb_y, nzt_y
	885	DO i = nxl_y, nxr_y
	886
	887	work(0:ny) = ar(0:ny,i,k)
	888	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	889
	890	DO j = 0, (ny+1)/2
	891	ar(j,i,k) = work(2*j)
	892	ENDDO
	893	DO j = 1, (ny+1)/2 - 1
	894	ar(ny+1-j,i,k) = work(2*j+1)
	895	ENDDO
	896
	897	ENDDO
	898	ENDDO
	899	!$OMP END PARALLEL
	900
	901	ELSE
	902
	903	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	904	!$OMP DO
	905	DO k = nzb_y, nzt_y
	906	DO i = nxl_y, nxr_y
	907
	908	DO j = 0, (ny+1)/2
	909	work(2*j) = ar(j,i,k)
	910	ENDDO
	911	DO j = 1, (ny+1)/2 - 1
	912	work(2*j+1) = ar(ny+1-j,i,k)
	913	ENDDO
	914	work(1) = 0.0
	915	work(ny+2) = 0.0
	916
	917	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	918	ar(0:ny,i,k) = work(0:ny)
	919
	920	ENDDO
	921	ENDDO
	922	!$OMP END PARALLEL
	923
	924	ENDIF
	925
	926	ELSEIF ( fft_method == 'system-specific' ) THEN
	927
	928	#if defined( __ibm ) && ! defined( __ibmy_special )
	929	IF ( forward_fft) THEN
	930
	931	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	932	!$OMP DO
	933	DO k = nzb_y, nzt_y
	934	DO i = nxl_y, nxr_y
	935
	936	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
	937	auy2, nau2 )
	938
	939	DO j = 0, (ny+1)/2
	940	ar(j,i,k) = work(2*j)
	941	ENDDO
	942	DO j = 1, (ny+1)/2 - 1
	943	ar(ny+1-j,i,k) = work(2*j+1)
	944	ENDDO
	945
	946	ENDDO
	947	ENDDO
	948	!$OMP END PARALLEL
	949
	950	ELSE
	951
	952	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	953	!$OMP DO
	954	DO k = nzb_y, nzt_y
	955	DO i = nxl_y, nxr_y
	956
	957	DO j = 0, (ny+1)/2
	958	work(2*j) = ar(j,i,k)
	959	ENDDO
	960	DO j = 1, (ny+1)/2 - 1
	961	work(2*j+1) = ar(ny+1-j,i,k)
	962	ENDDO
	963	work(1) = 0.0
	964	work(ny+2) = 0.0
	965
	966	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
	967	auy4, nau2 )
	968
	969	DO j = 0, ny
	970	ar(j,i,k) = work(j)
	971	ENDDO
	972
	973	ENDDO
	974	ENDDO
	975	!$OMP END PARALLEL
	976
	977	ENDIF
	978	#elif defined( __nec )
	979	IF ( forward_fft ) THEN
	980
	981	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	982	!$OMP DO
	983	DO k = nzb_y, nzt_y
	984	DO i = nxl_y, nxr_y
	985
	986	work(0:ny) = ar(0:ny,i,k)
	987
	988	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
	989
	990	DO j = 0, (ny+1)/2
	991	ar(j,i,k) = work(2*j)
	992	ENDDO
	993	DO j = 1, (ny+1)/2 - 1
	994	ar(ny+1-j,i,k) = work(2*j+1)
	995	ENDDO
	996
	997	ENDDO
	998	ENDDO
	999	!$END OMP PARALLEL
	1000
	1001	ELSE
	1002
	1003	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1004	!$OMP DO
	1005	DO k = nzb_y, nzt_y
	1006	DO i = nxl_y, nxr_y
	1007
	1008	DO j = 0, (ny+1)/2
	1009	work(2*j) = ar(j,i,k)
	1010	ENDDO
	1011	DO j = 1, (ny+1)/2 - 1
	1012	work(2*j+1) = ar(ny+1-j,i,k)
	1013	ENDDO
	1014	work(1) = 0.0
	1015	work(ny+2) = 0.0
	1016
	1017	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
	1018
	1019	ar(0:ny,i,k) = work(0:ny)
	1020
	1021	ENDDO
	1022	ENDDO
	1023	!$OMP END PARALLEL
	1024
	1025	ENDIF
	1026	#elif defined( __cuda_fft )
	1027
	1028	IF ( forward_fft ) THEN
	1029
[1111]	1030	!$acc data present( ar )
	1031	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
[1106]	1032
[1111]	1033	!$acc kernels
	1034	!$acc loop
[1106]	1035	DO k = nzb_y, nzt_y
	1036	DO i = nxl_y, nxr_y
	1037
[1111]	1038	!$acc loop vector( 32 )
[1106]	1039	DO j = 0, (ny+1)/2
[1111]	1040	ar(j,i,k) = REAL( ar_tmp(j,i,k) ) * dny
[1106]	1041	ENDDO
	1042
[1111]	1043	!$acc loop vector( 32 )
[1106]	1044	DO j = 1, (ny+1)/2 - 1
[1111]	1045	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
[1106]	1046	ENDDO
	1047
	1048	ENDDO
	1049	ENDDO
[1111]	1050	!$acc end kernels
	1051	!$acc end data
[1106]	1052
	1053	ELSE
	1054
[1111]	1055	!$acc data present( ar )
	1056	!$acc kernels
	1057	!$acc loop
[1106]	1058	DO k = nzb_y, nzt_y
	1059	DO i = nxl_y, nxr_y
	1060
[1111]	1061	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
[1106]	1062
[1111]	1063	!$acc loop vector( 32 )
[1106]	1064	DO j = 1, (ny+1)/2 - 1
[1111]	1065	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
[1106]	1066	ENDDO
[1111]	1067	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
[1106]	1068
	1069	ENDDO
	1070	ENDDO
[1111]	1071	!$acc end kernels
[1106]	1072
[1111]	1073	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
	1074	!$acc end data
[1106]	1075
	1076	ENDIF
	1077
	1078	#else
	1079	message_string = 'no system-specific fft-call available'
	1080	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
	1081	#endif
	1082
	1083	ELSE
	1084
	1085	message_string = 'fft method "' // TRIM( fft_method) // &
	1086	'" not available'
	1087	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
	1088
	1089	ENDIF
	1090
	1091	END SUBROUTINE fft_y
	1092
	1093	SUBROUTINE fft_y_1d( ar, direction )
	1094
	1095	!----------------------------------------------------------------------!
	1096	! fft_y_1d !
	1097	! !
	1098	! Fourier-transformation along y-direction !
	1099	! Version for 1D-decomposition !
	1100	! !
	1101	! fft_y uses internal algorithms (Singleton or Temperton) or !
	1102	! system-specific routines, if they are available !
	1103	!----------------------------------------------------------------------!
	1104
	1105	IMPLICIT NONE
	1106
	1107	CHARACTER (LEN=*) :: direction
[1]	1108	INTEGER :: j, jshape(1)
	1109
[1106]	1110	LOGICAL :: forward_fft
	1111
[1]	1112	REAL, DIMENSION(0:ny) :: ar
	1113	REAL, DIMENSION(0:ny+2) :: work
	1114	REAL, DIMENSION(ny+2) :: work1
	1115	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	1116	#if defined( __ibm )
	1117	REAL, DIMENSION(nau2) :: auy2, auy4
	1118	#elif defined( __nec )
	1119	REAL, DIMENSION(6*(ny+1)) :: work2
	1120	#endif
	1121
[1106]	1122	IF ( direction == 'forward' ) THEN
	1123	forward_fft = .TRUE.
	1124	ELSE
	1125	forward_fft = .FALSE.
	1126	ENDIF
	1127
[1]	1128	IF ( fft_method == 'singleton-algorithm' ) THEN
	1129
	1130	!
	1131	!-- Performing the fft with singleton's software works on every system,
	1132	!-- since it is part of the model
	1133	ALLOCATE( cwork(0:ny) )
	1134
[1106]	1135	IF ( forward_fft ) THEN
[1]	1136
	1137	DO j = 0, ny
	1138	cwork(j) = CMPLX( ar(j) )
	1139	ENDDO
	1140
	1141	jshape = SHAPE( cwork )
	1142	CALL FFTN( cwork, jshape )
	1143
	1144	DO j = 0, (ny+1)/2
	1145	ar(j) = REAL( cwork(j) )
	1146	ENDDO
	1147	DO j = 1, (ny+1)/2 - 1
	1148	ar(ny+1-j) = -AIMAG( cwork(j) )
	1149	ENDDO
	1150
	1151	ELSE
	1152
	1153	cwork(0) = CMPLX( ar(0), 0.0 )
	1154	DO j = 1, (ny+1)/2 - 1
	1155	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
	1156	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
	1157	ENDDO
	1158	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
	1159
	1160	jshape = SHAPE( cwork )
	1161	CALL FFTN( cwork, jshape, inv = .TRUE. )
	1162
	1163	DO j = 0, ny
	1164	ar(j) = REAL( cwork(j) )
	1165	ENDDO
	1166
	1167	ENDIF
	1168
	1169	DEALLOCATE( cwork )
	1170
	1171	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	1172
	1173	!
	1174	!-- Performing the fft with Temperton's software works on every system,
	1175	!-- since it is part of the model
[1106]	1176	IF ( forward_fft ) THEN
[1]	1177
	1178	work(0:ny) = ar
	1179	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	1180
	1181	DO j = 0, (ny+1)/2
	1182	ar(j) = work(2*j)
	1183	ENDDO
	1184	DO j = 1, (ny+1)/2 - 1
	1185	ar(ny+1-j) = work(2*j+1)
	1186	ENDDO
	1187
	1188	ELSE
	1189
	1190	DO j = 0, (ny+1)/2
	1191	work(2*j) = ar(j)
	1192	ENDDO
	1193	DO j = 1, (ny+1)/2 - 1
	1194	work(2*j+1) = ar(ny+1-j)
	1195	ENDDO
	1196	work(1) = 0.0
	1197	work(ny+2) = 0.0
	1198
	1199	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	1200	ar = work(0:ny)
	1201
	1202	ENDIF
	1203
	1204	ELSEIF ( fft_method == 'system-specific' ) THEN
	1205
	1206	#if defined( __ibm ) && ! defined( __ibmy_special )
[1106]	1207	IF ( forward_fft ) THEN
[1]	1208
[1106]	1209	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
[1]	1210	auy2, nau2 )
	1211
	1212	DO j = 0, (ny+1)/2
	1213	ar(j) = work(2*j)
	1214	ENDDO
	1215	DO j = 1, (ny+1)/2 - 1
	1216	ar(ny+1-j) = work(2*j+1)
	1217	ENDDO
	1218
	1219	ELSE
	1220
	1221	DO j = 0, (ny+1)/2
	1222	work(2*j) = ar(j)
	1223	ENDDO
	1224	DO j = 1, (ny+1)/2 - 1
	1225	work(2*j+1) = ar(ny+1-j)
	1226	ENDDO
	1227	work(1) = 0.0
	1228	work(ny+2) = 0.0
	1229
[1106]	1230	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]	1231	auy4, nau2 )
	1232
	1233	DO j = 0, ny
	1234	ar(j) = work(j)
	1235	ENDDO
	1236
	1237	ENDIF
	1238	#elif defined( __nec )
[1106]	1239	IF ( forward_fft ) THEN
[1]	1240
	1241	work(0:ny) = ar(0:ny)
	1242
[1106]	1243	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
[1]	1244
	1245	DO j = 0, (ny+1)/2
	1246	ar(j) = work(2*j)
	1247	ENDDO
	1248	DO j = 1, (ny+1)/2 - 1
	1249	ar(ny+1-j) = work(2*j+1)
	1250	ENDDO
	1251
	1252	ELSE
	1253
	1254	DO j = 0, (ny+1)/2
	1255	work(2*j) = ar(j)
	1256	ENDDO
	1257	DO j = 1, (ny+1)/2 - 1
	1258	work(2*j+1) = ar(ny+1-j)
	1259	ENDDO
	1260	work(1) = 0.0
	1261	work(ny+2) = 0.0
	1262
[1106]	1263	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
[1]	1264
	1265	ar(0:ny) = work(0:ny)
	1266
	1267	ENDIF
	1268	#else
[254]	1269	message_string = 'no system-specific fft-call available'
[1106]	1270	CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 )
[254]	1271
[1]	1272	#endif
	1273
	1274	ELSE
	1275
[274]	1276	message_string = 'fft method "' // TRIM( fft_method) // &
	1277	'" not available'
[1106]	1278	CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1279
	1280	ENDIF
	1281
[1106]	1282	END SUBROUTINE fft_y_1d
[1]	1283
	1284	SUBROUTINE fft_x_m( ar, direction )
	1285
	1286	!----------------------------------------------------------------------!
	1287	! fft_x_m !
	1288	! !
	1289	! Fourier-transformation along x-direction !
	1290	! Version for 1d domain decomposition !
	1291	! using multiple 1D FFT from Math Keisan on NEC !
	1292	! or Temperton-algorithm !
	1293	! (no singleton-algorithm on NEC because it does not vectorize) !
	1294	! !
	1295	!----------------------------------------------------------------------!
	1296
	1297	IMPLICIT NONE
	1298
	1299	CHARACTER (LEN=*) :: direction
[1092]	1300	INTEGER :: i, k, siza
[1]	1301
	1302	REAL, DIMENSION(0:nx,nz) :: ar
	1303	REAL, DIMENSION(0:nx+3,nz+1) :: ai
	1304	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
	1305	#if defined( __nec )
[1092]	1306	INTEGER :: sizw
[1]	1307	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
	1308	#endif
	1309
	1310	IF ( fft_method == 'temperton-algorithm' ) THEN
	1311
	1312	siza = SIZE( ai, 1 )
	1313
	1314	IF ( direction == 'forward') THEN
	1315
	1316	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	1317	ai(nx+1:,:) = 0.0
	1318
	1319	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
	1320
	1321	DO k = 1, nz
	1322	DO i = 0, (nx+1)/2
	1323	ar(i,k) = ai(2*i,k)
	1324	ENDDO
	1325	DO i = 1, (nx+1)/2 - 1
	1326	ar(nx+1-i,k) = ai(2*i+1,k)
	1327	ENDDO
	1328	ENDDO
	1329
	1330	ELSE
	1331
	1332	DO k = 1, nz
	1333	DO i = 0, (nx+1)/2
	1334	ai(2*i,k) = ar(i,k)
	1335	ENDDO
	1336	DO i = 1, (nx+1)/2 - 1
	1337	ai(2*i+1,k) = ar(nx+1-i,k)
	1338	ENDDO
	1339	ai(1,k) = 0.0
	1340	ai(nx+2,k) = 0.0
	1341	ENDDO
	1342
	1343	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
	1344
	1345	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	1346
	1347	ENDIF
	1348
	1349	ELSEIF ( fft_method == 'system-specific' ) THEN
	1350
	1351	#if defined( __nec )
	1352	siza = SIZE( ai, 1 )
	1353	sizw = SIZE( work, 1 )
	1354
	1355	IF ( direction == 'forward') THEN
	1356
	1357	!
	1358	!-- Tables are initialized once more. This call should not be
	1359	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1360	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]	1361	trig_xf, work1, 0 )
	1362
	1363	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	1364	IF ( nz1 > nz ) THEN
	1365	ai(:,nz1) = 0.0
	1366	ENDIF
	1367
[1106]	1368	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
[1]	1369	trig_xf, work1, 0 )
	1370
	1371	DO k = 1, nz
	1372	DO i = 0, (nx+1)/2
	1373	ar(i,k) = REAL( work(i+1,k) )
	1374	ENDDO
	1375	DO i = 1, (nx+1)/2 - 1
	1376	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
	1377	ENDDO
	1378	ENDDO
	1379
	1380	ELSE
	1381
	1382	!
	1383	!-- Tables are initialized once more. This call should not be
	1384	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1385	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]	1386	trig_xb, work1, 0 )
	1387
	1388	IF ( nz1 > nz ) THEN
	1389	work(:,nz1) = 0.0
	1390	ENDIF
	1391	DO k = 1, nz
	1392	work(1,k) = CMPLX( ar(0,k), 0.0 )
	1393	DO i = 1, (nx+1)/2 - 1
	1394	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
	1395	ENDDO
	1396	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
	1397	ENDDO
	1398
[1106]	1399	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
[1]	1400	trig_xb, work1, 0 )
	1401
	1402	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	1403
	1404	ENDIF
	1405
	1406	#else
[254]	1407	message_string = 'no system-specific fft-call available'
	1408	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	1409	#endif
	1410
	1411	ELSE
	1412
[274]	1413	message_string = 'fft method "' // TRIM( fft_method) // &
	1414	'" not available'
[254]	1415	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1416
	1417	ENDIF
	1418
	1419	END SUBROUTINE fft_x_m
	1420
	1421	SUBROUTINE fft_y_m( ar, ny1, direction )
	1422
	1423	!----------------------------------------------------------------------!
	1424	! fft_y_m !
	1425	! !
	1426	! Fourier-transformation along y-direction !
	1427	! Version for 1d domain decomposition !
	1428	! using multiple 1D FFT from Math Keisan on NEC !
	1429	! or Temperton-algorithm !
	1430	! (no singleton-algorithm on NEC because it does not vectorize) !
	1431	! !
	1432	!----------------------------------------------------------------------!
	1433
	1434	IMPLICIT NONE
	1435
	1436	CHARACTER (LEN=*) :: direction
[1092]	1437	INTEGER :: j, k, ny1, siza
[1]	1438
	1439	REAL, DIMENSION(0:ny1,nz) :: ar
	1440	REAL, DIMENSION(0:ny+3,nz+1) :: ai
	1441	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
	1442	#if defined( __nec )
[1092]	1443	INTEGER :: sizw
[1]	1444	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
	1445	#endif
	1446
	1447	IF ( fft_method == 'temperton-algorithm' ) THEN
	1448
	1449	siza = SIZE( ai, 1 )
	1450
	1451	IF ( direction == 'forward') THEN
	1452
	1453	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	1454	ai(ny+1:,:) = 0.0
	1455
	1456	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
	1457
	1458	DO k = 1, nz
	1459	DO j = 0, (ny+1)/2
	1460	ar(j,k) = ai(2*j,k)
	1461	ENDDO
	1462	DO j = 1, (ny+1)/2 - 1
	1463	ar(ny+1-j,k) = ai(2*j+1,k)
	1464	ENDDO
	1465	ENDDO
	1466
	1467	ELSE
	1468
	1469	DO k = 1, nz
	1470	DO j = 0, (ny+1)/2
	1471	ai(2*j,k) = ar(j,k)
	1472	ENDDO
	1473	DO j = 1, (ny+1)/2 - 1
	1474	ai(2*j+1,k) = ar(ny+1-j,k)
	1475	ENDDO
	1476	ai(1,k) = 0.0
	1477	ai(ny+2,k) = 0.0
	1478	ENDDO
	1479
	1480	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
	1481
	1482	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	1483
	1484	ENDIF
	1485
	1486	ELSEIF ( fft_method == 'system-specific' ) THEN
	1487
	1488	#if defined( __nec )
	1489	siza = SIZE( ai, 1 )
	1490	sizw = SIZE( work, 1 )
	1491
	1492	IF ( direction == 'forward') THEN
	1493
	1494	!
	1495	!-- Tables are initialized once more. This call should not be
	1496	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1497	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]	1498	trig_yf, work1, 0 )
	1499
	1500	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	1501	IF ( nz1 > nz ) THEN
	1502	ai(:,nz1) = 0.0
	1503	ENDIF
	1504
[1106]	1505	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
[1]	1506	trig_yf, work1, 0 )
	1507
	1508	DO k = 1, nz
	1509	DO j = 0, (ny+1)/2
	1510	ar(j,k) = REAL( work(j+1,k) )
	1511	ENDDO
	1512	DO j = 1, (ny+1)/2 - 1
	1513	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
	1514	ENDDO
	1515	ENDDO
	1516
	1517	ELSE
	1518
	1519	!
	1520	!-- Tables are initialized once more. This call should not be
	1521	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1522	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]	1523	trig_yb, work1, 0 )
	1524
	1525	IF ( nz1 > nz ) THEN
	1526	work(:,nz1) = 0.0
	1527	ENDIF
	1528	DO k = 1, nz
	1529	work(1,k) = CMPLX( ar(0,k), 0.0 )
	1530	DO j = 1, (ny+1)/2 - 1
	1531	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
	1532	ENDDO
	1533	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
	1534	ENDDO
	1535
[1106]	1536	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
[1]	1537	trig_yb, work1, 0 )
	1538
	1539	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	1540
	1541	ENDIF
	1542
	1543	#else
[254]	1544	message_string = 'no system-specific fft-call available'
	1545	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	1546	#endif
	1547
	1548	ELSE
[254]	1549
[274]	1550	message_string = 'fft method "' // TRIM( fft_method) // &
	1551	'" not available'
[254]	1552	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1553
	1554	ENDIF
	1555
	1556	END SUBROUTINE fft_y_m
	1557
[1106]	1558
[1]	1559	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |