Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |