Home

Context Navigation

fft_xy.f90 @ 1849

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

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

Context Navigation

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

Download in other formats: