Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |