Home

Context Navigation

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

Last change on this file since 1503 was 1483, checked in by raasch, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 58.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |