Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |