Home

Context Navigation

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

Last change on this file since 1762 was 1750, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |