Home

Context Navigation

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

Last change on this file since 1397 was 1393, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 58.2 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	!
[1310]	17	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1322]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[1375]	22	!
[1393]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: fft_xy.f90 1393 2014-05-06 09:14:41Z raasch $
	27	!
[1393]	28	! 1392 2014-05-06 09:10:05Z raasch
	29	! bugfix: KIND attribute added to CMPLX functions
	30	!
[1375]	31	! 1374 2014-04-25 12:55:07Z raasch
	32	! bugfixes: missing variables added to ONLY list, dpk renamed dp
	33	!
[1373]	34	! 1372 2014-04-24 06:29:32Z raasch
	35	! openMP-bugfix for fftw: some arrays defined as threadprivate
	36	!
[1354]	37	! 1353 2014-04-08 15:21:23Z heinze
	38	! REAL constants provided with KIND-attribute
	39	!
[1343]	40	! 1342 2014-03-26 17:04:47Z kanani
	41	! REAL constants defined as wp-kind
	42	!
[1323]	43	! 1322 2014-03-20 16:38:49Z raasch
	44	! REAL functions provided with KIND-attribute
	45	!
[1321]	46	! 1320 2014-03-20 08:40:49Z raasch
[1320]	47	! ONLY-attribute added to USE-statements,
	48	! kind-parameters added to all INTEGER and REAL declaration statements,
	49	! kinds are defined in new module kinds,
	50	! old module precision_kind is removed,
	51	! revision history before 2012 removed,
	52	! comment fields (!:) to be used for variable explanations added to
	53	! all variable declaration statements
[1]	54	!
[1305]	55	! 1304 2014-03-12 10:29:42Z raasch
	56	! openmp bugfix: work1 used in Temperton algorithm must be private
	57	!
[1258]	58	! 1257 2013-11-08 15:18:40Z raasch
	59	! openacc loop and loop vector clauses removed, declare create moved after
	60	! the FORTRAN declaration statement
	61	!
[1220]	62	! 1219 2013-08-30 09:33:18Z heinze
	63	! bugfix: use own branch for fftw
	64	!
[1217]	65	! 1216 2013-08-26 09:31:42Z raasch
	66	! fft_x and fft_y modified for parallel / ovverlapping execution of fft and
	67	! transpositions,
	68	! fftw implemented for 1d-decomposition (fft_x_1d, fft_y_1d)
	69	!
[1211]	70	! 1210 2013-08-14 10:58:20Z raasch
	71	! fftw added
	72	!
[1167]	73	! 1166 2013-05-24 13:55:44Z raasch
	74	! C_DOUBLE/COMPLEX reset to dpk
	75	!
[1154]	76	! 1153 2013-05-10 14:33:08Z raasch
	77	! code adjustment of data types for CUDA fft required by PGI 12.3 / CUDA 5.0
	78	!
[1112]	79	! 1111 2013-03-08 23:54:10Z raasch
	80	! further openACC statements added, CUDA branch completely runs on GPU
	81	! bugfix: CUDA fft plans adjusted for domain decomposition (before they always
	82	! used total domain)
	83	!
[1107]	84	! 1106 2013-03-04 05:31:38Z raasch
	85	! CUDA fft added
	86	! array_kind renamed precision_kind, 3D- instead of 1D-loops in fft_x and fft_y
	87	! old fft_x, fft_y become fft_x_1d, fft_y_1d and are used for 1D-decomposition
	88	!
[1093]	89	! 1092 2013-02-02 11:24:22Z raasch
	90	! variable sizw declared for NEC case only
	91	!
[1037]	92	! 1036 2012-10-22 13:43:42Z raasch
	93	! code put under GPL (PALM 3.9)
	94	!
[1]	95	! Revision 1.1 2002/06/11 13:00:49 raasch
	96	! Initial revision
	97	!
	98	!
	99	! Description:
	100	! ------------
	101	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
	102	! Original version: Klaus Ketelsen (May 2002)
	103	!------------------------------------------------------------------------------!
	104
[1320]	105	USE control_parameters, &
	106	ONLY: fft_method, message_string
	107
	108	USE indices, &
	109	ONLY: nx, ny, nz
	110
[1153]	111	#if defined( __cuda_fft )
	112	USE ISO_C_BINDING
[1210]	113	#elif defined( __fftw )
	114	USE, INTRINSIC :: ISO_C_BINDING
[1153]	115	#endif
[1320]	116
	117	USE kinds
	118
	119	USE singleton, &
	120	ONLY: fftn
	121
[1]	122	USE temperton_fft
[1320]	123
	124	USE transpose_indices, &
[1374]	125	ONLY: nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
[1]	126
	127	IMPLICIT NONE
	128
	129	PRIVATE
[1106]	130	PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
[1]	131
[1320]	132	INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x !:
	133	INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE :: ifax_y !:
[1]	134
[1320]	135	LOGICAL, SAVE :: init_fft = .FALSE. !:
[1]	136
[1320]	137	REAL(wp), SAVE :: dnx !:
	138	REAL(wp), SAVE :: dny !:
	139	REAL(wp), SAVE :: sqr_dnx !:
	140	REAL(wp), SAVE :: sqr_dny !:
	141
	142	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x !:
	143	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trigs_y !:
[1]	144
	145	#if defined( __ibm )
[1320]	146	INTEGER(iwp), PARAMETER :: nau1 = 20000 !:
	147	INTEGER(iwp), PARAMETER :: nau2 = 22000 !:
[1]	148	!
	149	!-- The following working arrays contain tables and have to be "save" and
	150	!-- shared in OpenMP sense
[1320]	151	REAL(wp), DIMENSION(nau1), SAVE :: aux1 !:
	152	REAL(wp), DIMENSION(nau1), SAVE :: auy1 !:
	153	REAL(wp), DIMENSION(nau1), SAVE :: aux3 !:
	154	REAL(wp), DIMENSION(nau1), SAVE :: auy3 !:
	155
[1]	156	#elif defined( __nec )
[1320]	157	INTEGER(iwp), SAVE :: nz1 !:
	158
	159	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb !:
	160	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_xf !:
	161	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_yb !:
	162	REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: trig_yf !:
	163
[1106]	164	#elif defined( __cuda_fft )
[1320]	165	INTEGER(C_INT), SAVE :: plan_xf !:
	166	INTEGER(C_INT), SAVE :: plan_xi !:
	167	INTEGER(C_INT), SAVE :: plan_yf !:
	168	INTEGER(C_INT), SAVE :: plan_yi !:
	169
	170	INTEGER(iwp), SAVE :: total_points_x_transpo !:
	171	INTEGER(iwp), SAVE :: total_points_y_transpo !:
[1219]	172	#endif
	173
	174	#if defined( __fftw )
[1210]	175	INCLUDE 'fftw3.f03'
[1320]	176	INTEGER(KIND=C_INT) :: nx_c !:
	177	INTEGER(KIND=C_INT) :: ny_c !:
	178
[1372]	179	!$OMP THREADPRIVATE( x_out, y_out, x_in, y_in )
[1320]	180	COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE :: &
	181	x_out !:
	182	COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE :: &
	183	y_out !:
	184
	185	REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE :: &
	186	x_in !:
	187	REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE :: &
	188	y_in !:
	189
	190
[1210]	191	TYPE(C_PTR), SAVE :: plan_xf, plan_xi, plan_yf, plan_yi
[1]	192	#endif
	193
	194	!
	195	!-- Public interfaces
	196	INTERFACE fft_init
	197	MODULE PROCEDURE fft_init
	198	END INTERFACE fft_init
	199
	200	INTERFACE fft_x
	201	MODULE PROCEDURE fft_x
	202	END INTERFACE fft_x
	203
[1106]	204	INTERFACE fft_x_1d
	205	MODULE PROCEDURE fft_x_1d
	206	END INTERFACE fft_x_1d
	207
[1]	208	INTERFACE fft_y
	209	MODULE PROCEDURE fft_y
	210	END INTERFACE fft_y
	211
[1106]	212	INTERFACE fft_y_1d
	213	MODULE PROCEDURE fft_y_1d
	214	END INTERFACE fft_y_1d
	215
[1]	216	INTERFACE fft_x_m
	217	MODULE PROCEDURE fft_x_m
	218	END INTERFACE fft_x_m
	219
	220	INTERFACE fft_y_m
	221	MODULE PROCEDURE fft_y_m
	222	END INTERFACE fft_y_m
	223
	224	CONTAINS
	225
	226
	227	SUBROUTINE fft_init
	228
[1106]	229	USE cuda_fft_interfaces
	230
[1]	231	IMPLICIT NONE
	232
	233	!
	234	!-- The following temporary working arrays have to be on stack or private
	235	!-- in OpenMP sense
	236	#if defined( __ibm )
[1320]	237	REAL(wp), DIMENSION(0:nx+2) :: workx !:
	238	REAL(wp), DIMENSION(0:ny+2) :: worky !:
	239	REAL(wp), DIMENSION(nau2) :: aux2 !:
	240	REAL(wp), DIMENSION(nau2) :: auy2 !:
	241	REAL(wp), DIMENSION(nau2) :: aux4 !:
	242	REAL(wp), DIMENSION(nau2) :: auy4 !:
[1]	243	#elif defined( __nec )
[1320]	244	REAL(wp), DIMENSION(0:nx+3,nz+1) :: work_x !:
	245	REAL(wp), DIMENSION(0:ny+3,nz+1) :: work_y !:
	246	REAL(wp), DIMENSION(6*(nx+3),nz+1) :: workx !:
	247	REAL(wp), DIMENSION(6*(ny+3),nz+1) :: worky !:
[1]	248	#endif
	249
	250	!
	251	!-- Return, if already called
	252	IF ( init_fft ) THEN
	253	RETURN
	254	ELSE
	255	init_fft = .TRUE.
	256	ENDIF
	257
	258	IF ( fft_method == 'system-specific' ) THEN
	259
[1342]	260	dnx = 1.0_wp / ( nx + 1.0_wp )
	261	dny = 1.0_wp / ( ny + 1.0_wp )
[1106]	262	sqr_dnx = SQRT( dnx )
	263	sqr_dny = SQRT( dny )
[1]	264	#if defined( __ibm ) && ! defined( __ibmy_special )
	265	!
	266	!-- Initialize tables for fft along x
[1106]	267	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
[1]	268	aux2, nau2 )
[1106]	269	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]	270	aux4, nau2 )
	271	!
	272	!-- Initialize tables for fft along y
[1106]	273	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
[1]	274	auy2, nau2 )
[1106]	275	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]	276	auy4, nau2 )
	277	#elif defined( __nec )
[254]	278	message_string = 'fft method "' // TRIM( fft_method) // &
	279	'" currently does not work on NEC'
	280	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
[1]	281
[1320]	282	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
[1]	283	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
	284
[1342]	285	work_x = 0.0_wp
	286	work_y = 0.0_wp
[1]	287	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
	288	! when using the NEC ffts
	289
	290	!
	291	!-- Initialize tables for fft along x (non-vector and vector case (M))
[1106]	292	CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
	293	CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
[1320]	294	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
[1]	295	trig_xf, workx, 0 )
[1320]	296	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, &
[1]	297	trig_xb, workx, 0 )
	298	!
	299	!-- Initialize tables for fft along y (non-vector and vector case (M))
[1106]	300	CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
	301	CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
[1320]	302	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
[1]	303	trig_yf, worky, 0 )
[1320]	304	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, &
[1]	305	trig_yb, worky, 0 )
[1106]	306	#elif defined( __cuda_fft )
	307	total_points_x_transpo = (nx+1) * (nyn_x-nys_x+1) * (nzt_x-nzb_x+1)
	308	total_points_y_transpo = (ny+1) * (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1)
[1111]	309	CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
	310	CALL CUFFTPLAN1D( plan_xi, nx+1, CUFFT_Z2D, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
	311	CALL CUFFTPLAN1D( plan_yf, ny+1, CUFFT_D2Z, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
	312	CALL CUFFTPLAN1D( plan_yi, ny+1, CUFFT_Z2D, (nxr_y-nxl_y+1) * (nzt_y-nzb_y+1) )
[1]	313	#else
[254]	314	message_string = 'no system-specific fft-call available'
	315	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
[1]	316	#endif
	317	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	318	!
	319	!-- Temperton-algorithm
	320	!-- Initialize tables for fft along x and y
	321	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
	322
	323	CALL set99( trigs_x, ifax_x, nx+1 )
	324	CALL set99( trigs_y, ifax_y, ny+1 )
	325
[1210]	326	ELSEIF ( fft_method == 'fftw' ) THEN
	327	!
	328	!-- FFTW
	329	#if defined( __fftw )
	330	nx_c = nx+1
	331	ny_c = ny+1
[1372]	332	!$OMP PARALLEL
[1320]	333	ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2), &
[1210]	334	y_out(0:(ny+1)/2) )
[1372]	335	!$OMP END PARALLEL
[1210]	336	plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
	337	plan_xi = FFTW_PLAN_DFT_C2R_1D( nx_c, x_out, x_in, FFTW_ESTIMATE )
	338	plan_yf = FFTW_PLAN_DFT_R2C_1D( ny_c, y_in, y_out, FFTW_ESTIMATE )
	339	plan_yi = FFTW_PLAN_DFT_C2R_1D( ny_c, y_out, y_in, FFTW_ESTIMATE )
	340	#else
	341	message_string = 'preprocessor switch for fftw is missing'
	342	CALL message( 'fft_init', 'PA0080', 1, 2, 0, 6, 0 )
	343	#endif
	344
[1]	345	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
	346
	347	CONTINUE
	348
	349	ELSE
	350
[254]	351	message_string = 'fft method "' // TRIM( fft_method) // &
	352	'" not available'
	353	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
[1]	354	ENDIF
	355
	356	END SUBROUTINE fft_init
	357
	358
[1216]	359	SUBROUTINE fft_x( ar, direction, ar_2d )
[1]	360
	361	!----------------------------------------------------------------------!
	362	! fft_x !
	363	! !
	364	! Fourier-transformation along x-direction !
[1106]	365	! Version for 2D-decomposition !
[1]	366	! !
	367	! fft_x uses internal algorithms (Singleton or Temperton) or !
	368	! system-specific routines, if they are available !
	369	!----------------------------------------------------------------------!
	370
[1106]	371	USE cuda_fft_interfaces
[1153]	372	#if defined( __cuda_fft )
	373	USE ISO_C_BINDING
	374	#endif
[1106]	375
[1]	376	IMPLICIT NONE
	377
[1320]	378	CHARACTER (LEN=*) :: direction !:
	379
	380	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !:
[1106]	381
[1320]	382	INTEGER(iwp) :: i !:
	383	INTEGER(iwp) :: ishape(1) !:
	384	INTEGER(iwp) :: j !:
	385	INTEGER(iwp) :: k !:
[1106]	386
[1320]	387	LOGICAL :: forward_fft !:
	388
	389	REAL(wp), DIMENSION(0:nx+2) :: work !:
	390	REAL(wp), DIMENSION(nx+2) :: work1 !:
	391
[1106]	392	#if defined( __ibm )
[1320]	393	REAL(wp), DIMENSION(nau2) :: aux2 !:
	394	REAL(wp), DIMENSION(nau2) :: aux4 !:
[1106]	395	#elif defined( __nec )
[1320]	396	REAL(wp), DIMENSION(6*(nx+1)) :: work2 !:
[1106]	397	#elif defined( __cuda_fft )
[1374]	398	COMPLEX(dp), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) :: &
[1320]	399	ar_tmp !:
[1111]	400	!$acc declare create( ar_tmp )
[1106]	401	#endif
	402
[1320]	403	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL :: &
	404	ar_2d !:
	405	REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: &
	406	ar !:
	407
[1106]	408	IF ( direction == 'forward' ) THEN
	409	forward_fft = .TRUE.
	410	ELSE
	411	forward_fft = .FALSE.
	412	ENDIF
	413
	414	IF ( fft_method == 'singleton-algorithm' ) THEN
	415
	416	!
	417	!-- Performing the fft with singleton's software works on every system,
	418	!-- since it is part of the model
	419	ALLOCATE( cwork(0:nx) )
	420
	421	IF ( forward_fft ) then
	422
	423	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
	424	!$OMP DO
	425	DO k = nzb_x, nzt_x
	426	DO j = nys_x, nyn_x
	427
	428	DO i = 0, nx
[1392]	429	cwork(i) = CMPLX( ar(i,j,k), KIND=wp )
[1106]	430	ENDDO
	431
	432	ishape = SHAPE( cwork )
	433	CALL FFTN( cwork, ishape )
	434
	435	DO i = 0, (nx+1)/2
[1322]	436	ar(i,j,k) = REAL( cwork(i), KIND=wp )
[1106]	437	ENDDO
	438	DO i = 1, (nx+1)/2 - 1
	439	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
	440	ENDDO
	441
	442	ENDDO
	443	ENDDO
	444	!$OMP END PARALLEL
	445
	446	ELSE
	447
	448	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
	449	!$OMP DO
	450	DO k = nzb_x, nzt_x
	451	DO j = nys_x, nyn_x
	452
[1392]	453	cwork(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
[1106]	454	DO i = 1, (nx+1)/2 - 1
[1392]	455	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k), &
	456	KIND=wp )
	457	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), &
	458	KIND=wp )
[1106]	459	ENDDO
[1392]	460	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
[1106]	461
	462	ishape = SHAPE( cwork )
	463	CALL FFTN( cwork, ishape, inv = .TRUE. )
	464
	465	DO i = 0, nx
[1322]	466	ar(i,j,k) = REAL( cwork(i), KIND=wp )
[1106]	467	ENDDO
	468
	469	ENDDO
	470	ENDDO
	471	!$OMP END PARALLEL
	472
	473	ENDIF
	474
	475	DEALLOCATE( cwork )
	476
	477	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	478
	479	!
	480	!-- Performing the fft with Temperton's software works on every system,
	481	!-- since it is part of the model
	482	IF ( forward_fft ) THEN
	483
[1304]	484	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]	485	!$OMP DO
	486	DO k = nzb_x, nzt_x
	487	DO j = nys_x, nyn_x
	488
	489	work(0:nx) = ar(0:nx,j,k)
	490	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	491
	492	DO i = 0, (nx+1)/2
	493	ar(i,j,k) = work(2*i)
	494	ENDDO
	495	DO i = 1, (nx+1)/2 - 1
	496	ar(nx+1-i,j,k) = work(2*i+1)
	497	ENDDO
	498
	499	ENDDO
	500	ENDDO
	501	!$OMP END PARALLEL
	502
	503	ELSE
	504
[1304]	505	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]	506	!$OMP DO
	507	DO k = nzb_x, nzt_x
	508	DO j = nys_x, nyn_x
	509
	510	DO i = 0, (nx+1)/2
	511	work(2*i) = ar(i,j,k)
	512	ENDDO
	513	DO i = 1, (nx+1)/2 - 1
	514	work(2*i+1) = ar(nx+1-i,j,k)
	515	ENDDO
[1342]	516	work(1) = 0.0_wp
	517	work(nx+2) = 0.0_wp
[1106]	518
	519	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	520	ar(0:nx,j,k) = work(0:nx)
	521
	522	ENDDO
	523	ENDDO
	524	!$OMP END PARALLEL
	525
	526	ENDIF
	527
[1210]	528	ELSEIF ( fft_method == 'fftw' ) THEN
	529
	530	#if defined( __fftw )
	531	IF ( forward_fft ) THEN
	532
	533	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	534	!$OMP DO
	535	DO k = nzb_x, nzt_x
	536	DO j = nys_x, nyn_x
	537
	538	x_in(0:nx) = ar(0:nx,j,k)
	539	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
	540
[1216]	541	IF ( PRESENT( ar_2d ) ) THEN
[1210]	542
[1216]	543	DO i = 0, (nx+1)/2
[1322]	544	ar_2d(i,j) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
[1216]	545	ENDDO
	546	DO i = 1, (nx+1)/2 - 1
	547	ar_2d(nx+1-i,j) = AIMAG( x_out(i) ) / ( nx+1 )
	548	ENDDO
	549
	550	ELSE
	551
	552	DO i = 0, (nx+1)/2
[1322]	553	ar(i,j,k) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
[1216]	554	ENDDO
	555	DO i = 1, (nx+1)/2 - 1
	556	ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
	557	ENDDO
	558
	559	ENDIF
	560
[1210]	561	ENDDO
	562	ENDDO
	563	!$OMP END PARALLEL
	564
[1216]	565	ELSE
[1210]	566	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	567	!$OMP DO
	568	DO k = nzb_x, nzt_x
	569	DO j = nys_x, nyn_x
	570
[1216]	571	IF ( PRESENT( ar_2d ) ) THEN
[1210]	572
[1392]	573	x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp, KIND=wp )
[1216]	574	DO i = 1, (nx+1)/2 - 1
[1392]	575	x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j), &
	576	KIND=wp )
[1216]	577	ENDDO
[1392]	578	x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp, &
	579	KIND=wp )
[1216]	580
	581	ELSE
	582
[1392]	583	x_out(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
[1216]	584	DO i = 1, (nx+1)/2 - 1
[1392]	585	x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
[1216]	586	ENDDO
[1392]	587	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, &
	588	KIND=wp )
[1216]	589
	590	ENDIF
	591
[1210]	592	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
	593	ar(0:nx,j,k) = x_in(0:nx)
	594
	595	ENDDO
	596	ENDDO
	597	!$OMP END PARALLEL
	598
[1216]	599	ENDIF
[1210]	600	#endif
	601
[1106]	602	ELSEIF ( fft_method == 'system-specific' ) THEN
	603
	604	#if defined( __ibm ) && ! defined( __ibmy_special )
	605	IF ( forward_fft ) THEN
	606
	607	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	608	!$OMP DO
	609	DO k = nzb_x, nzt_x
	610	DO j = nys_x, nyn_x
	611
[1320]	612	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, &
	613	nau1, aux2, nau2 )
[1106]	614
	615	DO i = 0, (nx+1)/2
	616	ar(i,j,k) = work(2*i)
	617	ENDDO
	618	DO i = 1, (nx+1)/2 - 1
	619	ar(nx+1-i,j,k) = work(2*i+1)
	620	ENDDO
	621
	622	ENDDO
	623	ENDDO
	624	!$OMP END PARALLEL
	625
	626	ELSE
	627
	628	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	629	!$OMP DO
	630	DO k = nzb_x, nzt_x
	631	DO j = nys_x, nyn_x
	632
	633	DO i = 0, (nx+1)/2
	634	work(2*i) = ar(i,j,k)
	635	ENDDO
	636	DO i = 1, (nx+1)/2 - 1
	637	work(2*i+1) = ar(nx+1-i,j,k)
	638	ENDDO
[1342]	639	work(1) = 0.0_wp
	640	work(nx+2) = 0.0_wp
[1106]	641
[1320]	642	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, &
	643	aux3, nau1, aux4, nau2 )
[1106]	644
	645	DO i = 0, nx
	646	ar(i,j,k) = work(i)
	647	ENDDO
	648
	649	ENDDO
	650	ENDDO
	651	!$OMP END PARALLEL
	652
	653	ENDIF
	654
	655	#elif defined( __nec )
	656
	657	IF ( forward_fft ) THEN
	658
	659	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	660	!$OMP DO
	661	DO k = nzb_x, nzt_x
	662	DO j = nys_x, nyn_x
	663
	664	work(0:nx) = ar(0:nx,j,k)
	665
	666	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
	667
	668	DO i = 0, (nx+1)/2
	669	ar(i,j,k) = work(2*i)
	670	ENDDO
	671	DO i = 1, (nx+1)/2 - 1
	672	ar(nx+1-i,j,k) = work(2*i+1)
	673	ENDDO
	674
	675	ENDDO
	676	ENDDO
	677	!$END OMP PARALLEL
	678
	679	ELSE
	680
	681	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	682	!$OMP DO
	683	DO k = nzb_x, nzt_x
	684	DO j = nys_x, nyn_x
	685
	686	DO i = 0, (nx+1)/2
	687	work(2*i) = ar(i,j,k)
	688	ENDDO
	689	DO i = 1, (nx+1)/2 - 1
	690	work(2*i+1) = ar(nx+1-i,j,k)
	691	ENDDO
[1342]	692	work(1) = 0.0_wp
	693	work(nx+2) = 0.0_wp
[1106]	694
	695	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
	696
	697	ar(0:nx,j,k) = work(0:nx)
	698
	699	ENDDO
	700	ENDDO
	701	!$OMP END PARALLEL
	702
	703	ENDIF
	704
	705	#elif defined( __cuda_fft )
	706
	707	IF ( forward_fft ) THEN
	708
[1111]	709	!$acc data present( ar )
	710	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
[1106]	711
[1111]	712	!$acc kernels
[1106]	713	DO k = nzb_x, nzt_x
	714	DO j = nys_x, nyn_x
	715
	716	DO i = 0, (nx+1)/2
[1322]	717	ar(i,j,k) = REAL( ar_tmp(i,j,k), KIND=wp ) * dnx
[1106]	718	ENDDO
	719
	720	DO i = 1, (nx+1)/2 - 1
[1111]	721	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
[1106]	722	ENDDO
	723
	724	ENDDO
	725	ENDDO
[1111]	726	!$acc end kernels
	727	!$acc end data
[1106]	728
	729	ELSE
	730
[1111]	731	!$acc data present( ar )
	732	!$acc kernels
[1106]	733	DO k = nzb_x, nzt_x
	734	DO j = nys_x, nyn_x
	735
[1392]	736	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
[1106]	737
	738	DO i = 1, (nx+1)/2 - 1
[1392]	739	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), &
	740	KIND=wp )
[1106]	741	ENDDO
[1392]	742	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, &
	743	KIND=wp )
[1106]	744
	745	ENDDO
	746	ENDDO
[1111]	747	!$acc end kernels
[1106]	748
[1111]	749	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
	750	!$acc end data
[1106]	751
	752	ENDIF
	753
	754	#else
	755	message_string = 'no system-specific fft-call available'
	756	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
	757	#endif
	758
	759	ELSE
	760
	761	message_string = 'fft method "' // TRIM( fft_method) // &
	762	'" not available'
	763	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
	764
	765	ENDIF
	766
	767	END SUBROUTINE fft_x
	768
	769	SUBROUTINE fft_x_1d( ar, direction )
	770
	771	!----------------------------------------------------------------------!
	772	! fft_x_1d !
	773	! !
	774	! Fourier-transformation along x-direction !
	775	! Version for 1D-decomposition !
	776	! !
	777	! fft_x uses internal algorithms (Singleton or Temperton) or !
	778	! system-specific routines, if they are available !
	779	!----------------------------------------------------------------------!
	780
	781	IMPLICIT NONE
	782
[1320]	783	CHARACTER (LEN=*) :: direction !:
	784
	785	INTEGER(iwp) :: i !:
	786	INTEGER(iwp) :: ishape(1) !:
[1]	787
[1320]	788	LOGICAL :: forward_fft !:
[1106]	789
[1320]	790	REAL(wp), DIMENSION(0:nx) :: ar !:
	791	REAL(wp), DIMENSION(0:nx+2) :: work !:
	792	REAL(wp), DIMENSION(nx+2) :: work1 !:
	793
	794	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !:
	795
[1]	796	#if defined( __ibm )
[1320]	797	REAL(wp), DIMENSION(nau2) :: aux2 !:
	798	REAL(wp), DIMENSION(nau2) :: aux4 !:
[1]	799	#elif defined( __nec )
[1320]	800	REAL(wp), DIMENSION(6*(nx+1)) :: work2 !:
[1]	801	#endif
	802
[1106]	803	IF ( direction == 'forward' ) THEN
	804	forward_fft = .TRUE.
	805	ELSE
	806	forward_fft = .FALSE.
	807	ENDIF
	808
[1]	809	IF ( fft_method == 'singleton-algorithm' ) THEN
	810
	811	!
	812	!-- Performing the fft with singleton's software works on every system,
	813	!-- since it is part of the model
	814	ALLOCATE( cwork(0:nx) )
	815
[1106]	816	IF ( forward_fft ) then
[1]	817
	818	DO i = 0, nx
[1392]	819	cwork(i) = CMPLX( ar(i), KIND=wp )
[1]	820	ENDDO
	821	ishape = SHAPE( cwork )
	822	CALL FFTN( cwork, ishape )
	823	DO i = 0, (nx+1)/2
[1322]	824	ar(i) = REAL( cwork(i), KIND=wp )
[1]	825	ENDDO
	826	DO i = 1, (nx+1)/2 - 1
	827	ar(nx+1-i) = -AIMAG( cwork(i) )
	828	ENDDO
	829
	830	ELSE
	831
[1392]	832	cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1]	833	DO i = 1, (nx+1)/2 - 1
[1392]	834	cwork(i) = CMPLX( ar(i), -ar(nx+1-i), KIND=wp )
	835	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i), KIND=wp )
[1]	836	ENDDO
[1392]	837	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
[1]	838
	839	ishape = SHAPE( cwork )
	840	CALL FFTN( cwork, ishape, inv = .TRUE. )
	841
	842	DO i = 0, nx
[1322]	843	ar(i) = REAL( cwork(i), KIND=wp )
[1]	844	ENDDO
	845
	846	ENDIF
	847
	848	DEALLOCATE( cwork )
	849
	850	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	851
	852	!
	853	!-- Performing the fft with Temperton's software works on every system,
	854	!-- since it is part of the model
[1106]	855	IF ( forward_fft ) THEN
[1]	856
	857	work(0:nx) = ar
	858	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	859
	860	DO i = 0, (nx+1)/2
	861	ar(i) = work(2*i)
	862	ENDDO
	863	DO i = 1, (nx+1)/2 - 1
	864	ar(nx+1-i) = work(2*i+1)
	865	ENDDO
	866
	867	ELSE
	868
	869	DO i = 0, (nx+1)/2
	870	work(2*i) = ar(i)
	871	ENDDO
	872	DO i = 1, (nx+1)/2 - 1
	873	work(2*i+1) = ar(nx+1-i)
	874	ENDDO
[1342]	875	work(1) = 0.0_wp
	876	work(nx+2) = 0.0_wp
[1]	877
	878	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	879	ar = work(0:nx)
	880
	881	ENDIF
	882
[1216]	883	ELSEIF ( fft_method == 'fftw' ) THEN
	884
	885	#if defined( __fftw )
	886	IF ( forward_fft ) THEN
	887
	888	x_in(0:nx) = ar(0:nx)
	889	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
	890
	891	DO i = 0, (nx+1)/2
[1322]	892	ar(i) = REAL( x_out(i), KIND=wp ) / ( nx+1 )
[1216]	893	ENDDO
	894	DO i = 1, (nx+1)/2 - 1
	895	ar(nx+1-i) = AIMAG( x_out(i) ) / ( nx+1 )
	896	ENDDO
	897
	898	ELSE
	899
[1392]	900	x_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1216]	901	DO i = 1, (nx+1)/2 - 1
[1392]	902	x_out(i) = CMPLX( ar(i), ar(nx+1-i), KIND=wp )
[1216]	903	ENDDO
[1392]	904	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0_wp, KIND=wp )
[1216]	905
	906	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
	907	ar(0:nx) = x_in(0:nx)
	908
	909	ENDIF
	910	#endif
	911
[1]	912	ELSEIF ( fft_method == 'system-specific' ) THEN
	913
	914	#if defined( __ibm ) && ! defined( __ibmy_special )
[1106]	915	IF ( forward_fft ) THEN
[1]	916
[1320]	917	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
[1]	918	aux2, nau2 )
	919
	920	DO i = 0, (nx+1)/2
	921	ar(i) = work(2*i)
	922	ENDDO
	923	DO i = 1, (nx+1)/2 - 1
	924	ar(nx+1-i) = work(2*i+1)
	925	ENDDO
	926
	927	ELSE
	928
	929	DO i = 0, (nx+1)/2
	930	work(2*i) = ar(i)
	931	ENDDO
	932	DO i = 1, (nx+1)/2 - 1
	933	work(2*i+1) = ar(nx+1-i)
	934	ENDDO
[1342]	935	work(1) = 0.0_wp
	936	work(nx+2) = 0.0_wp
[1]	937
[1106]	938	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]	939	aux4, nau2 )
	940
	941	DO i = 0, nx
	942	ar(i) = work(i)
	943	ENDDO
	944
	945	ENDIF
	946	#elif defined( __nec )
[1106]	947	IF ( forward_fft ) THEN
[1]	948
	949	work(0:nx) = ar(0:nx)
	950
[1106]	951	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
	952
[1]	953	DO i = 0, (nx+1)/2
	954	ar(i) = work(2*i)
	955	ENDDO
	956	DO i = 1, (nx+1)/2 - 1
	957	ar(nx+1-i) = work(2*i+1)
	958	ENDDO
	959
	960	ELSE
	961
	962	DO i = 0, (nx+1)/2
	963	work(2*i) = ar(i)
	964	ENDDO
	965	DO i = 1, (nx+1)/2 - 1
	966	work(2*i+1) = ar(nx+1-i)
	967	ENDDO
[1342]	968	work(1) = 0.0_wp
	969	work(nx+2) = 0.0_wp
[1]	970
[1106]	971	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
[1]	972
	973	ar(0:nx) = work(0:nx)
	974
	975	ENDIF
	976	#else
[254]	977	message_string = 'no system-specific fft-call available'
[1106]	978	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
[1]	979	#endif
	980	ELSE
[274]	981	message_string = 'fft method "' // TRIM( fft_method) // &
	982	'" not available'
[1106]	983	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]	984
	985	ENDIF
	986
[1106]	987	END SUBROUTINE fft_x_1d
[1]	988
[1216]	989	SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
	990	nxr_y_l )
[1]	991
	992	!----------------------------------------------------------------------!
	993	! fft_y !
	994	! !
	995	! Fourier-transformation along y-direction !
[1106]	996	! Version for 2D-decomposition !
[1]	997	! !
	998	! fft_y uses internal algorithms (Singleton or Temperton) or !
	999	! system-specific routines, if they are available !
[1216]	1000	! !
	1001	! direction: 'forward' or 'backward' !
	1002	! ar, ar_tr: 3D data arrays !
	1003	! forward: ar: before ar_tr: after transformation !
	1004	! backward: ar_tr: before ar: after transfosition !
	1005	! !
	1006	! In case of non-overlapping transposition/transformation: !
	1007	! nxl_y_bound = nxl_y_l = nxl_y !
	1008	! nxr_y_bound = nxr_y_l = nxr_y !
	1009	! !
	1010	! In case of overlapping transposition/transformation !
	1011	! - nxl_y_bound and nxr_y_bound have the original values of !
	1012	! nxl_y, nxr_y. ar_tr is dimensioned using these values. !
	1013	! - nxl_y_l = nxr_y_r. ar is dimensioned with these values, so that !
	1014	! transformation is carried out for a 2D-plane only. !
[1]	1015	!----------------------------------------------------------------------!
	1016
[1106]	1017	USE cuda_fft_interfaces
[1153]	1018	#if defined( __cuda_fft )
	1019	USE ISO_C_BINDING
	1020	#endif
[1106]	1021
[1]	1022	IMPLICIT NONE
	1023
[1320]	1024	CHARACTER (LEN=*) :: direction !:
	1025
	1026	INTEGER(iwp) :: i !:
	1027	INTEGER(iwp) :: j !:
	1028	INTEGER(iwp) :: jshape(1) !:
	1029	INTEGER(iwp) :: k !:
	1030	INTEGER(iwp) :: nxl_y_bound !:
	1031	INTEGER(iwp) :: nxl_y_l !:
	1032	INTEGER(iwp) :: nxr_y_bound !:
	1033	INTEGER(iwp) :: nxr_y_l !:
[1106]	1034
[1320]	1035	LOGICAL :: forward_fft !:
[1106]	1036
[1320]	1037	REAL(wp), DIMENSION(0:ny+2) :: work !:
	1038	REAL(wp), DIMENSION(ny+2) :: work1 !:
	1039
	1040	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !:
	1041
[1106]	1042	#if defined( __ibm )
[1320]	1043	REAL(wp), DIMENSION(nau2) :: auy2 !:
	1044	REAL(wp), DIMENSION(nau2) :: auy4 !:
[1106]	1045	#elif defined( __nec )
[1320]	1046	REAL(wp), DIMENSION(6*(ny+1)) :: work2 !:
[1106]	1047	#elif defined( __cuda_fft )
[1374]	1048	COMPLEX(dp), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: &
[1320]	1049	ar_tmp !:
[1111]	1050	!$acc declare create( ar_tmp )
[1106]	1051	#endif
	1052
[1320]	1053	REAL(wp), DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y) :: &
	1054	ar !:
	1055	REAL(wp), DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) :: &
	1056	ar_tr !:
	1057
[1106]	1058	IF ( direction == 'forward' ) THEN
	1059	forward_fft = .TRUE.
	1060	ELSE
	1061	forward_fft = .FALSE.
	1062	ENDIF
	1063
	1064	IF ( fft_method == 'singleton-algorithm' ) THEN
	1065
	1066	!
	1067	!-- Performing the fft with singleton's software works on every system,
	1068	!-- since it is part of the model
	1069	ALLOCATE( cwork(0:ny) )
	1070
	1071	IF ( forward_fft ) then
	1072
	1073	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
	1074	!$OMP DO
	1075	DO k = nzb_y, nzt_y
[1216]	1076	DO i = nxl_y_l, nxr_y_l
[1106]	1077
	1078	DO j = 0, ny
[1392]	1079	cwork(j) = CMPLX( ar(j,i,k), KIND=wp )
[1106]	1080	ENDDO
	1081
	1082	jshape = SHAPE( cwork )
	1083	CALL FFTN( cwork, jshape )
	1084
	1085	DO j = 0, (ny+1)/2
[1322]	1086	ar_tr(j,i,k) = REAL( cwork(j), KIND=wp )
[1106]	1087	ENDDO
	1088	DO j = 1, (ny+1)/2 - 1
[1216]	1089	ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
[1106]	1090	ENDDO
	1091
	1092	ENDDO
	1093	ENDDO
	1094	!$OMP END PARALLEL
	1095
	1096	ELSE
	1097
	1098	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
	1099	!$OMP DO
	1100	DO k = nzb_y, nzt_y
[1216]	1101	DO i = nxl_y_l, nxr_y_l
[1106]	1102
[1392]	1103	cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
[1106]	1104	DO j = 1, (ny+1)/2 - 1
[1392]	1105	cwork(j) = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k), &
	1106	KIND=wp )
	1107	cwork(ny+1-j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k), &
	1108	KIND=wp )
[1106]	1109	ENDDO
[1392]	1110	cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp, &
	1111	KIND=wp )
[1106]	1112
	1113	jshape = SHAPE( cwork )
	1114	CALL FFTN( cwork, jshape, inv = .TRUE. )
	1115
	1116	DO j = 0, ny
[1322]	1117	ar(j,i,k) = REAL( cwork(j), KIND=wp )
[1106]	1118	ENDDO
	1119
	1120	ENDDO
	1121	ENDDO
	1122	!$OMP END PARALLEL
	1123
	1124	ENDIF
	1125
	1126	DEALLOCATE( cwork )
	1127
	1128	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	1129
	1130	!
	1131	!-- Performing the fft with Temperton's software works on every system,
	1132	!-- since it is part of the model
	1133	IF ( forward_fft ) THEN
	1134
[1304]	1135	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]	1136	!$OMP DO
	1137	DO k = nzb_y, nzt_y
[1216]	1138	DO i = nxl_y_l, nxr_y_l
[1106]	1139
	1140	work(0:ny) = ar(0:ny,i,k)
	1141	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	1142
	1143	DO j = 0, (ny+1)/2
[1216]	1144	ar_tr(j,i,k) = work(2*j)
[1106]	1145	ENDDO
	1146	DO j = 1, (ny+1)/2 - 1
[1216]	1147	ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]	1148	ENDDO
	1149
	1150	ENDDO
	1151	ENDDO
	1152	!$OMP END PARALLEL
	1153
	1154	ELSE
	1155
[1304]	1156	!$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
[1106]	1157	!$OMP DO
	1158	DO k = nzb_y, nzt_y
[1216]	1159	DO i = nxl_y_l, nxr_y_l
[1106]	1160
	1161	DO j = 0, (ny+1)/2
[1216]	1162	work(2*j) = ar_tr(j,i,k)
[1106]	1163	ENDDO
	1164	DO j = 1, (ny+1)/2 - 1
[1216]	1165	work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]	1166	ENDDO
[1342]	1167	work(1) = 0.0_wp
	1168	work(ny+2) = 0.0_wp
[1106]	1169
	1170	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	1171	ar(0:ny,i,k) = work(0:ny)
	1172
	1173	ENDDO
	1174	ENDDO
	1175	!$OMP END PARALLEL
	1176
	1177	ENDIF
	1178
[1210]	1179	ELSEIF ( fft_method == 'fftw' ) THEN
	1180
	1181	#if defined( __fftw )
	1182	IF ( forward_fft ) THEN
	1183
	1184	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1185	!$OMP DO
	1186	DO k = nzb_y, nzt_y
[1216]	1187	DO i = nxl_y_l, nxr_y_l
[1210]	1188
	1189	y_in(0:ny) = ar(0:ny,i,k)
	1190	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
	1191
	1192	DO j = 0, (ny+1)/2
[1322]	1193	ar_tr(j,i,k) = REAL( y_out(j), KIND=wp ) / (ny+1)
[1210]	1194	ENDDO
	1195	DO j = 1, (ny+1)/2 - 1
[1216]	1196	ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
[1210]	1197	ENDDO
	1198
	1199	ENDDO
	1200	ENDDO
	1201	!$OMP END PARALLEL
	1202
	1203	ELSE
	1204
	1205	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1206	!$OMP DO
	1207	DO k = nzb_y, nzt_y
[1216]	1208	DO i = nxl_y_l, nxr_y_l
[1210]	1209
[1392]	1210	y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
[1210]	1211	DO j = 1, (ny+1)/2 - 1
[1392]	1212	y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k, &
	1213	KIND=wp) )
[1210]	1214	ENDDO
[1392]	1215	y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp, &
	1216	KIND=wp )
[1210]	1217
	1218	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
	1219	ar(0:ny,i,k) = y_in(0:ny)
	1220
	1221	ENDDO
	1222	ENDDO
	1223	!$OMP END PARALLEL
	1224
	1225	ENDIF
	1226	#endif
	1227
[1106]	1228	ELSEIF ( fft_method == 'system-specific' ) THEN
	1229
	1230	#if defined( __ibm ) && ! defined( __ibmy_special )
	1231	IF ( forward_fft) THEN
	1232
	1233	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1234	!$OMP DO
	1235	DO k = nzb_y, nzt_y
[1216]	1236	DO i = nxl_y_l, nxr_y_l
[1106]	1237
[1320]	1238	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, &
	1239	nau1, auy2, nau2 )
[1106]	1240
	1241	DO j = 0, (ny+1)/2
[1216]	1242	ar_tr(j,i,k) = work(2*j)
[1106]	1243	ENDDO
	1244	DO j = 1, (ny+1)/2 - 1
[1216]	1245	ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]	1246	ENDDO
	1247
	1248	ENDDO
	1249	ENDDO
	1250	!$OMP END PARALLEL
	1251
	1252	ELSE
	1253
	1254	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1255	!$OMP DO
	1256	DO k = nzb_y, nzt_y
[1216]	1257	DO i = nxl_y_l, nxr_y_l
[1106]	1258
	1259	DO j = 0, (ny+1)/2
[1216]	1260	work(2*j) = ar_tr(j,i,k)
[1106]	1261	ENDDO
	1262	DO j = 1, (ny+1)/2 - 1
[1216]	1263	work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]	1264	ENDDO
[1342]	1265	work(1) = 0.0_wp
	1266	work(ny+2) = 0.0_wp
[1106]	1267
[1320]	1268	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, &
	1269	auy3, nau1, auy4, nau2 )
[1106]	1270
	1271	DO j = 0, ny
	1272	ar(j,i,k) = work(j)
	1273	ENDDO
	1274
	1275	ENDDO
	1276	ENDDO
	1277	!$OMP END PARALLEL
	1278
	1279	ENDIF
	1280	#elif defined( __nec )
	1281	IF ( forward_fft ) THEN
	1282
	1283	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1284	!$OMP DO
	1285	DO k = nzb_y, nzt_y
[1216]	1286	DO i = nxl_y_l, nxr_y_l
[1106]	1287
	1288	work(0:ny) = ar(0:ny,i,k)
	1289
	1290	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
	1291
	1292	DO j = 0, (ny+1)/2
[1216]	1293	ar_tr(j,i,k) = work(2*j)
[1106]	1294	ENDDO
	1295	DO j = 1, (ny+1)/2 - 1
[1216]	1296	ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]	1297	ENDDO
	1298
	1299	ENDDO
	1300	ENDDO
	1301	!$END OMP PARALLEL
	1302
	1303	ELSE
	1304
	1305	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1306	!$OMP DO
	1307	DO k = nzb_y, nzt_y
[1216]	1308	DO i = nxl_y_l, nxr_y_l
[1106]	1309
	1310	DO j = 0, (ny+1)/2
[1216]	1311	work(2*j) = ar_tr(j,i,k)
[1106]	1312	ENDDO
	1313	DO j = 1, (ny+1)/2 - 1
[1216]	1314	work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]	1315	ENDDO
[1342]	1316	work(1) = 0.0_wp
	1317	work(ny+2) = 0.0_wp
[1106]	1318
	1319	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
	1320
	1321	ar(0:ny,i,k) = work(0:ny)
	1322
	1323	ENDDO
	1324	ENDDO
	1325	!$OMP END PARALLEL
	1326
	1327	ENDIF
	1328	#elif defined( __cuda_fft )
	1329
	1330	IF ( forward_fft ) THEN
	1331
[1111]	1332	!$acc data present( ar )
	1333	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
[1106]	1334
[1111]	1335	!$acc kernels
[1106]	1336	DO k = nzb_y, nzt_y
	1337	DO i = nxl_y, nxr_y
	1338
	1339	DO j = 0, (ny+1)/2
[1322]	1340	ar(j,i,k) = REAL( ar_tmp(j,i,k), KIND=wp ) * dny
[1106]	1341	ENDDO
	1342
	1343	DO j = 1, (ny+1)/2 - 1
[1111]	1344	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
[1106]	1345	ENDDO
	1346
	1347	ENDDO
	1348	ENDDO
[1111]	1349	!$acc end kernels
	1350	!$acc end data
[1106]	1351
	1352	ELSE
	1353
[1111]	1354	!$acc data present( ar )
	1355	!$acc kernels
[1106]	1356	DO k = nzb_y, nzt_y
	1357	DO i = nxl_y, nxr_y
	1358
[1392]	1359	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0_wp, KIND=wp )
[1106]	1360
	1361	DO j = 1, (ny+1)/2 - 1
[1392]	1362	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k), &
	1363	KIND=wp )
[1106]	1364	ENDDO
[1392]	1365	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp, &
	1366	KIND=wp )
[1106]	1367
	1368	ENDDO
	1369	ENDDO
[1111]	1370	!$acc end kernels
[1106]	1371
[1111]	1372	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
	1373	!$acc end data
[1106]	1374
	1375	ENDIF
	1376
	1377	#else
	1378	message_string = 'no system-specific fft-call available'
	1379	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
	1380	#endif
	1381
	1382	ELSE
	1383
	1384	message_string = 'fft method "' // TRIM( fft_method) // &
	1385	'" not available'
	1386	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
	1387
	1388	ENDIF
	1389
	1390	END SUBROUTINE fft_y
	1391
	1392	SUBROUTINE fft_y_1d( ar, direction )
	1393
	1394	!----------------------------------------------------------------------!
	1395	! fft_y_1d !
	1396	! !
	1397	! Fourier-transformation along y-direction !
	1398	! Version for 1D-decomposition !
	1399	! !
	1400	! fft_y uses internal algorithms (Singleton or Temperton) or !
	1401	! system-specific routines, if they are available !
	1402	!----------------------------------------------------------------------!
	1403
	1404	IMPLICIT NONE
	1405
	1406	CHARACTER (LEN=*) :: direction
[1320]	1407
	1408	INTEGER(iwp) :: j !:
	1409	INTEGER(iwp) :: jshape(1) !:
[1]	1410
[1320]	1411	LOGICAL :: forward_fft !:
[1106]	1412
[1320]	1413	REAL(wp), DIMENSION(0:ny) :: ar !:
	1414	REAL(wp), DIMENSION(0:ny+2) :: work !:
	1415	REAL(wp), DIMENSION(ny+2) :: work1 !:
	1416
	1417	COMPLEX(wp), DIMENSION(:), ALLOCATABLE :: cwork !:
	1418
[1]	1419	#if defined( __ibm )
[1320]	1420	REAL(wp), DIMENSION(nau2) :: auy2 !:
	1421	REAL(wp), DIMENSION(nau2) :: auy4 !:
[1]	1422	#elif defined( __nec )
[1320]	1423	REAL(wp), DIMENSION(6*(ny+1)) :: work2 !:
[1]	1424	#endif
	1425
[1106]	1426	IF ( direction == 'forward' ) THEN
	1427	forward_fft = .TRUE.
	1428	ELSE
	1429	forward_fft = .FALSE.
	1430	ENDIF
	1431
[1]	1432	IF ( fft_method == 'singleton-algorithm' ) THEN
	1433
	1434	!
	1435	!-- Performing the fft with singleton's software works on every system,
	1436	!-- since it is part of the model
	1437	ALLOCATE( cwork(0:ny) )
	1438
[1106]	1439	IF ( forward_fft ) THEN
[1]	1440
	1441	DO j = 0, ny
[1392]	1442	cwork(j) = CMPLX( ar(j), KIND=wp )
[1]	1443	ENDDO
	1444
	1445	jshape = SHAPE( cwork )
	1446	CALL FFTN( cwork, jshape )
	1447
	1448	DO j = 0, (ny+1)/2
[1322]	1449	ar(j) = REAL( cwork(j), KIND=wp )
[1]	1450	ENDDO
	1451	DO j = 1, (ny+1)/2 - 1
	1452	ar(ny+1-j) = -AIMAG( cwork(j) )
	1453	ENDDO
	1454
	1455	ELSE
	1456
[1392]	1457	cwork(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1]	1458	DO j = 1, (ny+1)/2 - 1
[1392]	1459	cwork(j) = CMPLX( ar(j), -ar(ny+1-j), KIND=wp )
	1460	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j), KIND=wp )
[1]	1461	ENDDO
[1392]	1462	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
[1]	1463
	1464	jshape = SHAPE( cwork )
	1465	CALL FFTN( cwork, jshape, inv = .TRUE. )
	1466
	1467	DO j = 0, ny
[1322]	1468	ar(j) = REAL( cwork(j), KIND=wp )
[1]	1469	ENDDO
	1470
	1471	ENDIF
	1472
	1473	DEALLOCATE( cwork )
	1474
	1475	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	1476
	1477	!
	1478	!-- Performing the fft with Temperton's software works on every system,
	1479	!-- since it is part of the model
[1106]	1480	IF ( forward_fft ) THEN
[1]	1481
	1482	work(0:ny) = ar
	1483	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	1484
	1485	DO j = 0, (ny+1)/2
	1486	ar(j) = work(2*j)
	1487	ENDDO
	1488	DO j = 1, (ny+1)/2 - 1
	1489	ar(ny+1-j) = work(2*j+1)
	1490	ENDDO
	1491
	1492	ELSE
	1493
	1494	DO j = 0, (ny+1)/2
	1495	work(2*j) = ar(j)
	1496	ENDDO
	1497	DO j = 1, (ny+1)/2 - 1
	1498	work(2*j+1) = ar(ny+1-j)
	1499	ENDDO
[1342]	1500	work(1) = 0.0_wp
	1501	work(ny+2) = 0.0_wp
[1]	1502
	1503	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	1504	ar = work(0:ny)
	1505
	1506	ENDIF
	1507
[1216]	1508	ELSEIF ( fft_method == 'fftw' ) THEN
	1509
	1510	#if defined( __fftw )
	1511	IF ( forward_fft ) THEN
	1512
	1513	y_in(0:ny) = ar(0:ny)
	1514	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
	1515
	1516	DO j = 0, (ny+1)/2
[1322]	1517	ar(j) = REAL( y_out(j), KIND=wp ) / (ny+1)
[1216]	1518	ENDDO
	1519	DO j = 1, (ny+1)/2 - 1
	1520	ar(ny+1-j) = AIMAG( y_out(j) ) / (ny+1)
	1521	ENDDO
	1522
	1523	ELSE
	1524
[1392]	1525	y_out(0) = CMPLX( ar(0), 0.0_wp, KIND=wp )
[1216]	1526	DO j = 1, (ny+1)/2 - 1
[1392]	1527	y_out(j) = CMPLX( ar(j), ar(ny+1-j), KIND=wp )
[1216]	1528	ENDDO
[1392]	1529	y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0_wp, KIND=wp )
[1216]	1530
	1531	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
	1532	ar(0:ny) = y_in(0:ny)
	1533
	1534	ENDIF
	1535	#endif
	1536
[1]	1537	ELSEIF ( fft_method == 'system-specific' ) THEN
	1538
	1539	#if defined( __ibm ) && ! defined( __ibmy_special )
[1106]	1540	IF ( forward_fft ) THEN
[1]	1541
[1320]	1542	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
[1]	1543	auy2, nau2 )
	1544
	1545	DO j = 0, (ny+1)/2
	1546	ar(j) = work(2*j)
	1547	ENDDO
	1548	DO j = 1, (ny+1)/2 - 1
	1549	ar(ny+1-j) = work(2*j+1)
	1550	ENDDO
	1551
	1552	ELSE
	1553
	1554	DO j = 0, (ny+1)/2
	1555	work(2*j) = ar(j)
	1556	ENDDO
	1557	DO j = 1, (ny+1)/2 - 1
	1558	work(2*j+1) = ar(ny+1-j)
	1559	ENDDO
[1342]	1560	work(1) = 0.0_wp
	1561	work(ny+2) = 0.0_wp
[1]	1562
[1320]	1563	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, &
	1564	nau1, auy4, nau2 )
[1]	1565
	1566	DO j = 0, ny
	1567	ar(j) = work(j)
	1568	ENDDO
	1569
	1570	ENDIF
	1571	#elif defined( __nec )
[1106]	1572	IF ( forward_fft ) THEN
[1]	1573
	1574	work(0:ny) = ar(0:ny)
	1575
[1106]	1576	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
[1]	1577
	1578	DO j = 0, (ny+1)/2
	1579	ar(j) = work(2*j)
	1580	ENDDO
	1581	DO j = 1, (ny+1)/2 - 1
	1582	ar(ny+1-j) = work(2*j+1)
	1583	ENDDO
	1584
	1585	ELSE
	1586
	1587	DO j = 0, (ny+1)/2
	1588	work(2*j) = ar(j)
	1589	ENDDO
	1590	DO j = 1, (ny+1)/2 - 1
	1591	work(2*j+1) = ar(ny+1-j)
	1592	ENDDO
[1342]	1593	work(1) = 0.0_wp
	1594	work(ny+2) = 0.0_wp
[1]	1595
[1106]	1596	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
[1]	1597
	1598	ar(0:ny) = work(0:ny)
	1599
	1600	ENDIF
	1601	#else
[254]	1602	message_string = 'no system-specific fft-call available'
[1106]	1603	CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 )
[254]	1604
[1]	1605	#endif
	1606
	1607	ELSE
	1608
[274]	1609	message_string = 'fft method "' // TRIM( fft_method) // &
	1610	'" not available'
[1106]	1611	CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1612
	1613	ENDIF
	1614
[1106]	1615	END SUBROUTINE fft_y_1d
[1]	1616
	1617	SUBROUTINE fft_x_m( ar, direction )
	1618
	1619	!----------------------------------------------------------------------!
	1620	! fft_x_m !
	1621	! !
	1622	! Fourier-transformation along x-direction !
	1623	! Version for 1d domain decomposition !
	1624	! using multiple 1D FFT from Math Keisan on NEC !
	1625	! or Temperton-algorithm !
	1626	! (no singleton-algorithm on NEC because it does not vectorize) !
	1627	! !
	1628	!----------------------------------------------------------------------!
	1629
	1630	IMPLICIT NONE
	1631
[1320]	1632	CHARACTER (LEN=*) :: direction !:
	1633
	1634	INTEGER(iwp) :: i !:
	1635	INTEGER(iwp) :: k !:
	1636	INTEGER(iwp) :: siza !:
[1]	1637
[1320]	1638	REAL(wp), DIMENSION(0:nx,nz) :: ar !:
	1639	REAL(wp), DIMENSION(0:nx+3,nz+1) :: ai !:
	1640	REAL(wp), DIMENSION(6*(nx+4),nz+1) :: work1 !:
	1641
[1]	1642	#if defined( __nec )
[1320]	1643	INTEGER(iwp) :: sizw !:
	1644
	1645	COMPLEX(wp), DIMENSION((nx+4)/2+1,nz+1) :: work !:
[1]	1646	#endif
	1647
	1648	IF ( fft_method == 'temperton-algorithm' ) THEN
	1649
	1650	siza = SIZE( ai, 1 )
	1651
	1652	IF ( direction == 'forward') THEN
	1653
	1654	ai(0:nx,1:nz) = ar(0:nx,1:nz)
[1342]	1655	ai(nx+1:,:) = 0.0_wp
[1]	1656
	1657	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
	1658
	1659	DO k = 1, nz
	1660	DO i = 0, (nx+1)/2
	1661	ar(i,k) = ai(2*i,k)
	1662	ENDDO
	1663	DO i = 1, (nx+1)/2 - 1
	1664	ar(nx+1-i,k) = ai(2*i+1,k)
	1665	ENDDO
	1666	ENDDO
	1667
	1668	ELSE
	1669
	1670	DO k = 1, nz
	1671	DO i = 0, (nx+1)/2
	1672	ai(2*i,k) = ar(i,k)
	1673	ENDDO
	1674	DO i = 1, (nx+1)/2 - 1
	1675	ai(2*i+1,k) = ar(nx+1-i,k)
	1676	ENDDO
[1342]	1677	ai(1,k) = 0.0_wp
	1678	ai(nx+2,k) = 0.0_wp
[1]	1679	ENDDO
	1680
	1681	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
	1682
	1683	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	1684
	1685	ENDIF
	1686
	1687	ELSEIF ( fft_method == 'system-specific' ) THEN
	1688
	1689	#if defined( __nec )
	1690	siza = SIZE( ai, 1 )
	1691	sizw = SIZE( work, 1 )
	1692
	1693	IF ( direction == 'forward') THEN
	1694
	1695	!
	1696	!-- Tables are initialized once more. This call should not be
	1697	!-- necessary, but otherwise program aborts in asymmetric case
[1320]	1698	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]	1699	trig_xf, work1, 0 )
	1700
	1701	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	1702	IF ( nz1 > nz ) THEN
[1342]	1703	ai(:,nz1) = 0.0_wp
[1]	1704	ENDIF
	1705
[1320]	1706	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
[1]	1707	trig_xf, work1, 0 )
	1708
	1709	DO k = 1, nz
	1710	DO i = 0, (nx+1)/2
[1322]	1711	ar(i,k) = REAL( work(i+1,k), KIND=wp )
[1]	1712	ENDDO
	1713	DO i = 1, (nx+1)/2 - 1
	1714	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
	1715	ENDDO
	1716	ENDDO
	1717
	1718	ELSE
	1719
	1720	!
	1721	!-- Tables are initialized once more. This call should not be
	1722	!-- necessary, but otherwise program aborts in asymmetric case
[1320]	1723	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]	1724	trig_xb, work1, 0 )
	1725
	1726	IF ( nz1 > nz ) THEN
[1342]	1727	work(:,nz1) = 0.0_wp
[1]	1728	ENDIF
	1729	DO k = 1, nz
[1392]	1730	work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
[1]	1731	DO i = 1, (nx+1)/2 - 1
[1392]	1732	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k), KIND=wp )
[1]	1733	ENDDO
[1392]	1734	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0_wp, KIND=wp )
[1]	1735	ENDDO
	1736
[1106]	1737	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
[1]	1738	trig_xb, work1, 0 )
	1739
	1740	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	1741
	1742	ENDIF
	1743
	1744	#else
[254]	1745	message_string = 'no system-specific fft-call available'
	1746	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	1747	#endif
	1748
	1749	ELSE
	1750
[274]	1751	message_string = 'fft method "' // TRIM( fft_method) // &
	1752	'" not available'
[254]	1753	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1754
	1755	ENDIF
	1756
	1757	END SUBROUTINE fft_x_m
	1758
	1759	SUBROUTINE fft_y_m( ar, ny1, direction )
	1760
	1761	!----------------------------------------------------------------------!
	1762	! fft_y_m !
	1763	! !
	1764	! Fourier-transformation along y-direction !
	1765	! Version for 1d domain decomposition !
	1766	! using multiple 1D FFT from Math Keisan on NEC !
	1767	! or Temperton-algorithm !
	1768	! (no singleton-algorithm on NEC because it does not vectorize) !
	1769	! !
	1770	!----------------------------------------------------------------------!
	1771
	1772	IMPLICIT NONE
	1773
[1320]	1774	CHARACTER (LEN=*) :: direction !:
	1775
	1776	INTEGER(iwp) :: j !:
	1777	INTEGER(iwp) :: k !:
	1778	INTEGER(iwp) :: ny1 !:
	1779	INTEGER(iwp) :: siza !:
[1]	1780
[1320]	1781	REAL(wp), DIMENSION(0:ny1,nz) :: ar !:
	1782	REAL(wp), DIMENSION(0:ny+3,nz+1) :: ai !:
	1783	REAL(wp), DIMENSION(6*(ny+4),nz+1) :: work1 !:
	1784
[1]	1785	#if defined( __nec )
[1320]	1786	INTEGER(iwp) :: sizw !:
	1787
	1788	COMPLEX(wp), DIMENSION((ny+4)/2+1,nz+1) :: work !:
[1]	1789	#endif
	1790
	1791	IF ( fft_method == 'temperton-algorithm' ) THEN
	1792
	1793	siza = SIZE( ai, 1 )
	1794
	1795	IF ( direction == 'forward') THEN
	1796
	1797	ai(0:ny,1:nz) = ar(0:ny,1:nz)
[1342]	1798	ai(ny+1:,:) = 0.0_wp
[1]	1799
	1800	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
	1801
	1802	DO k = 1, nz
	1803	DO j = 0, (ny+1)/2
	1804	ar(j,k) = ai(2*j,k)
	1805	ENDDO
	1806	DO j = 1, (ny+1)/2 - 1
	1807	ar(ny+1-j,k) = ai(2*j+1,k)
	1808	ENDDO
	1809	ENDDO
	1810
	1811	ELSE
	1812
	1813	DO k = 1, nz
	1814	DO j = 0, (ny+1)/2
	1815	ai(2*j,k) = ar(j,k)
	1816	ENDDO
	1817	DO j = 1, (ny+1)/2 - 1
	1818	ai(2*j+1,k) = ar(ny+1-j,k)
	1819	ENDDO
[1342]	1820	ai(1,k) = 0.0_wp
	1821	ai(ny+2,k) = 0.0_wp
[1]	1822	ENDDO
	1823
	1824	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
	1825
	1826	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	1827
	1828	ENDIF
	1829
	1830	ELSEIF ( fft_method == 'system-specific' ) THEN
	1831
	1832	#if defined( __nec )
	1833	siza = SIZE( ai, 1 )
	1834	sizw = SIZE( work, 1 )
	1835
	1836	IF ( direction == 'forward') THEN
	1837
	1838	!
	1839	!-- Tables are initialized once more. This call should not be
	1840	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1841	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]	1842	trig_yf, work1, 0 )
	1843
	1844	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	1845	IF ( nz1 > nz ) THEN
[1342]	1846	ai(:,nz1) = 0.0_wp
[1]	1847	ENDIF
	1848
[1106]	1849	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
[1]	1850	trig_yf, work1, 0 )
	1851
	1852	DO k = 1, nz
	1853	DO j = 0, (ny+1)/2
[1322]	1854	ar(j,k) = REAL( work(j+1,k), KIND=wp )
[1]	1855	ENDDO
	1856	DO j = 1, (ny+1)/2 - 1
	1857	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
	1858	ENDDO
	1859	ENDDO
	1860
	1861	ELSE
	1862
	1863	!
	1864	!-- Tables are initialized once more. This call should not be
	1865	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1866	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]	1867	trig_yb, work1, 0 )
	1868
	1869	IF ( nz1 > nz ) THEN
[1342]	1870	work(:,nz1) = 0.0_wp
[1]	1871	ENDIF
	1872	DO k = 1, nz
[1392]	1873	work(1,k) = CMPLX( ar(0,k), 0.0_wp, KIND=wp )
[1]	1874	DO j = 1, (ny+1)/2 - 1
[1392]	1875	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k), KIND=wp )
[1]	1876	ENDDO
[1392]	1877	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0_wp, KIND=wp )
[1]	1878	ENDDO
	1879
[1106]	1880	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
[1]	1881	trig_yb, work1, 0 )
	1882
	1883	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	1884
	1885	ENDIF
	1886
	1887	#else
[254]	1888	message_string = 'no system-specific fft-call available'
	1889	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	1890	#endif
	1891
	1892	ELSE
[254]	1893
[274]	1894	message_string = 'fft method "' // TRIM( fft_method) // &
	1895	'" not available'
[254]	1896	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1897
	1898	ENDIF
	1899
	1900	END SUBROUTINE fft_y_m
	1901
[1106]	1902
[1]	1903	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |