Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |