Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |