Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |