Home

Context Navigation

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

Last change on this file since 1214 was 1211, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 49.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |