Home

Context Navigation

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

Last change on this file since 1859 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |