Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |