Home

Context Navigation

source: palm/trunk/SOURCE/fft_xy_mod.f90 @ 2106

Last change on this file since 2106 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |