Home

Context Navigation

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

Last change on this file since 1216 was 1216, checked in by raasch, 11 years ago
overlapping execution of fft and transpositions (MPI_ALLTOALL), but real overlapping is not activated so far, fftw implemented for 1D-decomposition resorting of arrays moved to separate routines resort_for_... bugfix in mbuild concerning Makefile_check
Property svn:keywords set to `Id`
File size: 53.5 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 1216 2013-08-26 09:31:42Z raasch $
[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
[1216]	288	SUBROUTINE fft_x( ar, direction, ar_2d )
[1]	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
[1216]	323	REAL, DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL :: ar_2d
[1106]	324	REAL, DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) :: ar
	325
	326	IF ( direction == 'forward' ) THEN
	327	forward_fft = .TRUE.
	328	ELSE
	329	forward_fft = .FALSE.
	330	ENDIF
	331
	332	IF ( fft_method == 'singleton-algorithm' ) THEN
	333
	334	!
	335	!-- Performing the fft with singleton's software works on every system,
	336	!-- since it is part of the model
	337	ALLOCATE( cwork(0:nx) )
	338
	339	IF ( forward_fft ) then
	340
	341	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
	342	!$OMP DO
	343	DO k = nzb_x, nzt_x
	344	DO j = nys_x, nyn_x
	345
	346	DO i = 0, nx
	347	cwork(i) = CMPLX( ar(i,j,k) )
	348	ENDDO
	349
	350	ishape = SHAPE( cwork )
	351	CALL FFTN( cwork, ishape )
	352
	353	DO i = 0, (nx+1)/2
	354	ar(i,j,k) = REAL( cwork(i) )
	355	ENDDO
	356	DO i = 1, (nx+1)/2 - 1
	357	ar(nx+1-i,j,k) = -AIMAG( cwork(i) )
	358	ENDDO
	359
	360	ENDDO
	361	ENDDO
	362	!$OMP END PARALLEL
	363
	364	ELSE
	365
	366	!$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
	367	!$OMP DO
	368	DO k = nzb_x, nzt_x
	369	DO j = nys_x, nyn_x
	370
	371	cwork(0) = CMPLX( ar(0,j,k), 0.0 )
	372	DO i = 1, (nx+1)/2 - 1
	373	cwork(i) = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k) )
	374	cwork(nx+1-i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
	375	ENDDO
	376	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
	377
	378	ishape = SHAPE( cwork )
	379	CALL FFTN( cwork, ishape, inv = .TRUE. )
	380
	381	DO i = 0, nx
	382	ar(i,j,k) = REAL( cwork(i) )
	383	ENDDO
	384
	385	ENDDO
	386	ENDDO
	387	!$OMP END PARALLEL
	388
	389	ENDIF
	390
	391	DEALLOCATE( cwork )
	392
	393	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	394
	395	!
	396	!-- Performing the fft with Temperton's software works on every system,
	397	!-- since it is part of the model
	398	IF ( forward_fft ) THEN
	399
	400	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	401	!$OMP DO
	402	DO k = nzb_x, nzt_x
	403	DO j = nys_x, nyn_x
	404
	405	work(0:nx) = ar(0:nx,j,k)
	406	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	407
	408	DO i = 0, (nx+1)/2
	409	ar(i,j,k) = work(2*i)
	410	ENDDO
	411	DO i = 1, (nx+1)/2 - 1
	412	ar(nx+1-i,j,k) = work(2*i+1)
	413	ENDDO
	414
	415	ENDDO
	416	ENDDO
	417	!$OMP END PARALLEL
	418
	419	ELSE
	420
	421	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	422	!$OMP DO
	423	DO k = nzb_x, nzt_x
	424	DO j = nys_x, nyn_x
	425
	426	DO i = 0, (nx+1)/2
	427	work(2*i) = ar(i,j,k)
	428	ENDDO
	429	DO i = 1, (nx+1)/2 - 1
	430	work(2*i+1) = ar(nx+1-i,j,k)
	431	ENDDO
	432	work(1) = 0.0
	433	work(nx+2) = 0.0
	434
	435	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	436	ar(0:nx,j,k) = work(0:nx)
	437
	438	ENDDO
	439	ENDDO
	440	!$OMP END PARALLEL
	441
	442	ENDIF
	443
[1210]	444	ELSEIF ( fft_method == 'fftw' ) THEN
	445
	446	#if defined( __fftw )
	447	IF ( forward_fft ) THEN
	448
	449	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	450	!$OMP DO
	451	DO k = nzb_x, nzt_x
	452	DO j = nys_x, nyn_x
	453
	454	x_in(0:nx) = ar(0:nx,j,k)
	455	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
	456
[1216]	457	IF ( PRESENT( ar_2d ) ) THEN
[1210]	458
[1216]	459	DO i = 0, (nx+1)/2
	460	ar_2d(i,j) = REAL( x_out(i) ) / ( nx+1 )
	461	ENDDO
	462	DO i = 1, (nx+1)/2 - 1
	463	ar_2d(nx+1-i,j) = AIMAG( x_out(i) ) / ( nx+1 )
	464	ENDDO
	465
	466	ELSE
	467
	468	DO i = 0, (nx+1)/2
	469	ar(i,j,k) = REAL( x_out(i) ) / ( nx+1 )
	470	ENDDO
	471	DO i = 1, (nx+1)/2 - 1
	472	ar(nx+1-i,j,k) = AIMAG( x_out(i) ) / ( nx+1 )
	473	ENDDO
	474
	475	ENDIF
	476
[1210]	477	ENDDO
	478	ENDDO
	479	!$OMP END PARALLEL
	480
[1216]	481	ELSE
[1210]	482	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	483	!$OMP DO
	484	DO k = nzb_x, nzt_x
	485	DO j = nys_x, nyn_x
	486
[1216]	487	IF ( PRESENT( ar_2d ) ) THEN
[1210]	488
[1216]	489	x_out(0) = CMPLX( ar_2d(0,j), 0.0 )
	490	DO i = 1, (nx+1)/2 - 1
	491	x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j) )
	492	ENDDO
	493	x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0 )
	494
	495	ELSE
	496
	497	x_out(0) = CMPLX( ar(0,j,k), 0.0 )
	498	DO i = 1, (nx+1)/2 - 1
	499	x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
	500	ENDDO
	501	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
	502
	503	ENDIF
	504
[1210]	505	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
	506	ar(0:nx,j,k) = x_in(0:nx)
	507
	508	ENDDO
	509	ENDDO
	510	!$OMP END PARALLEL
	511
[1216]	512	ENDIF
[1210]	513	#endif
	514
[1106]	515	ELSEIF ( fft_method == 'system-specific' ) THEN
	516
	517	#if defined( __ibm ) && ! defined( __ibmy_special )
	518	IF ( forward_fft ) THEN
	519
	520	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	521	!$OMP DO
	522	DO k = nzb_x, nzt_x
	523	DO j = nys_x, nyn_x
	524
	525	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
	526	aux2, nau2 )
	527
	528	DO i = 0, (nx+1)/2
	529	ar(i,j,k) = work(2*i)
	530	ENDDO
	531	DO i = 1, (nx+1)/2 - 1
	532	ar(nx+1-i,j,k) = work(2*i+1)
	533	ENDDO
	534
	535	ENDDO
	536	ENDDO
	537	!$OMP END PARALLEL
	538
	539	ELSE
	540
	541	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	542	!$OMP DO
	543	DO k = nzb_x, nzt_x
	544	DO j = nys_x, nyn_x
	545
	546	DO i = 0, (nx+1)/2
	547	work(2*i) = ar(i,j,k)
	548	ENDDO
	549	DO i = 1, (nx+1)/2 - 1
	550	work(2*i+1) = ar(nx+1-i,j,k)
	551	ENDDO
	552	work(1) = 0.0
	553	work(nx+2) = 0.0
	554
	555	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
	556	aux4, nau2 )
	557
	558	DO i = 0, nx
	559	ar(i,j,k) = work(i)
	560	ENDDO
	561
	562	ENDDO
	563	ENDDO
	564	!$OMP END PARALLEL
	565
	566	ENDIF
	567
	568	#elif defined( __nec )
	569
	570	IF ( forward_fft ) THEN
	571
	572	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	573	!$OMP DO
	574	DO k = nzb_x, nzt_x
	575	DO j = nys_x, nyn_x
	576
	577	work(0:nx) = ar(0:nx,j,k)
	578
	579	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
	580
	581	DO i = 0, (nx+1)/2
	582	ar(i,j,k) = work(2*i)
	583	ENDDO
	584	DO i = 1, (nx+1)/2 - 1
	585	ar(nx+1-i,j,k) = work(2*i+1)
	586	ENDDO
	587
	588	ENDDO
	589	ENDDO
	590	!$END OMP PARALLEL
	591
	592	ELSE
	593
	594	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	595	!$OMP DO
	596	DO k = nzb_x, nzt_x
	597	DO j = nys_x, nyn_x
	598
	599	DO i = 0, (nx+1)/2
	600	work(2*i) = ar(i,j,k)
	601	ENDDO
	602	DO i = 1, (nx+1)/2 - 1
	603	work(2*i+1) = ar(nx+1-i,j,k)
	604	ENDDO
	605	work(1) = 0.0
	606	work(nx+2) = 0.0
	607
	608	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
	609
	610	ar(0:nx,j,k) = work(0:nx)
	611
	612	ENDDO
	613	ENDDO
	614	!$OMP END PARALLEL
	615
	616	ENDIF
	617
	618	#elif defined( __cuda_fft )
	619
	620	IF ( forward_fft ) THEN
	621
[1111]	622	!$acc data present( ar )
	623	CALL CUFFTEXECD2Z( plan_xf, ar, ar_tmp )
[1106]	624
[1111]	625	!$acc kernels
	626	!$acc loop
[1106]	627	DO k = nzb_x, nzt_x
	628	DO j = nys_x, nyn_x
	629
[1111]	630	!$acc loop vector( 32 )
[1106]	631	DO i = 0, (nx+1)/2
[1111]	632	ar(i,j,k) = REAL( ar_tmp(i,j,k) ) * dnx
[1106]	633	ENDDO
	634
[1111]	635	!$acc loop vector( 32 )
[1106]	636	DO i = 1, (nx+1)/2 - 1
[1111]	637	ar(nx+1-i,j,k) = AIMAG( ar_tmp(i,j,k) ) * dnx
[1106]	638	ENDDO
	639
	640	ENDDO
	641	ENDDO
[1111]	642	!$acc end kernels
	643	!$acc end data
[1106]	644
	645	ELSE
	646
[1111]	647	!$acc data present( ar )
	648	!$acc kernels
	649	!$acc loop
[1106]	650	DO k = nzb_x, nzt_x
	651	DO j = nys_x, nyn_x
	652
[1111]	653	ar_tmp(0,j,k) = CMPLX( ar(0,j,k), 0.0 )
[1106]	654
[1111]	655	!$acc loop vector( 32 )
[1106]	656	DO i = 1, (nx+1)/2 - 1
[1111]	657	ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k) )
[1106]	658	ENDDO
[1111]	659	ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0 )
[1106]	660
	661	ENDDO
	662	ENDDO
[1111]	663	!$acc end kernels
[1106]	664
[1111]	665	CALL CUFFTEXECZ2D( plan_xi, ar_tmp, ar )
	666	!$acc end data
[1106]	667
	668	ENDIF
	669
	670	#else
	671	message_string = 'no system-specific fft-call available'
	672	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
	673	#endif
	674
	675	ELSE
	676
	677	message_string = 'fft method "' // TRIM( fft_method) // &
	678	'" not available'
	679	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
	680
	681	ENDIF
	682
	683	END SUBROUTINE fft_x
	684
	685	SUBROUTINE fft_x_1d( ar, direction )
	686
	687	!----------------------------------------------------------------------!
	688	! fft_x_1d !
	689	! !
	690	! Fourier-transformation along x-direction !
	691	! Version for 1D-decomposition !
	692	! !
	693	! fft_x uses internal algorithms (Singleton or Temperton) or !
	694	! system-specific routines, if they are available !
	695	!----------------------------------------------------------------------!
	696
	697	IMPLICIT NONE
	698
	699	CHARACTER (LEN=*) :: direction
[1]	700	INTEGER :: i, ishape(1)
	701
[1106]	702	LOGICAL :: forward_fft
	703
[1]	704	REAL, DIMENSION(0:nx) :: ar
	705	REAL, DIMENSION(0:nx+2) :: work
	706	REAL, DIMENSION(nx+2) :: work1
	707	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	708	#if defined( __ibm )
	709	REAL, DIMENSION(nau2) :: aux2, aux4
	710	#elif defined( __nec )
	711	REAL, DIMENSION(6*(nx+1)) :: work2
	712	#endif
	713
[1106]	714	IF ( direction == 'forward' ) THEN
	715	forward_fft = .TRUE.
	716	ELSE
	717	forward_fft = .FALSE.
	718	ENDIF
	719
[1]	720	IF ( fft_method == 'singleton-algorithm' ) THEN
	721
	722	!
	723	!-- Performing the fft with singleton's software works on every system,
	724	!-- since it is part of the model
	725	ALLOCATE( cwork(0:nx) )
	726
[1106]	727	IF ( forward_fft ) then
[1]	728
	729	DO i = 0, nx
	730	cwork(i) = CMPLX( ar(i) )
	731	ENDDO
	732	ishape = SHAPE( cwork )
	733	CALL FFTN( cwork, ishape )
	734	DO i = 0, (nx+1)/2
	735	ar(i) = REAL( cwork(i) )
	736	ENDDO
	737	DO i = 1, (nx+1)/2 - 1
	738	ar(nx+1-i) = -AIMAG( cwork(i) )
	739	ENDDO
	740
	741	ELSE
	742
	743	cwork(0) = CMPLX( ar(0), 0.0 )
	744	DO i = 1, (nx+1)/2 - 1
	745	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
	746	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
	747	ENDDO
	748	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
	749
	750	ishape = SHAPE( cwork )
	751	CALL FFTN( cwork, ishape, inv = .TRUE. )
	752
	753	DO i = 0, nx
	754	ar(i) = REAL( cwork(i) )
	755	ENDDO
	756
	757	ENDIF
	758
	759	DEALLOCATE( cwork )
	760
	761	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	762
	763	!
	764	!-- Performing the fft with Temperton's software works on every system,
	765	!-- since it is part of the model
[1106]	766	IF ( forward_fft ) THEN
[1]	767
	768	work(0:nx) = ar
	769	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	770
	771	DO i = 0, (nx+1)/2
	772	ar(i) = work(2*i)
	773	ENDDO
	774	DO i = 1, (nx+1)/2 - 1
	775	ar(nx+1-i) = work(2*i+1)
	776	ENDDO
	777
	778	ELSE
	779
	780	DO i = 0, (nx+1)/2
	781	work(2*i) = ar(i)
	782	ENDDO
	783	DO i = 1, (nx+1)/2 - 1
	784	work(2*i+1) = ar(nx+1-i)
	785	ENDDO
	786	work(1) = 0.0
	787	work(nx+2) = 0.0
	788
	789	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	790	ar = work(0:nx)
	791
	792	ENDIF
	793
[1216]	794	ELSEIF ( fft_method == 'fftw' ) THEN
	795
	796	#if defined( __fftw )
	797	IF ( forward_fft ) THEN
	798
	799	x_in(0:nx) = ar(0:nx)
	800	CALL FFTW_EXECUTE_DFT_R2C( plan_xf, x_in, x_out )
	801
	802	DO i = 0, (nx+1)/2
	803	ar(i) = REAL( x_out(i) ) / ( nx+1 )
	804	ENDDO
	805	DO i = 1, (nx+1)/2 - 1
	806	ar(nx+1-i) = AIMAG( x_out(i) ) / ( nx+1 )
	807	ENDDO
	808
	809	ELSE
	810
	811	x_out(0) = CMPLX( ar(0), 0.0 )
	812	DO i = 1, (nx+1)/2 - 1
	813	x_out(i) = CMPLX( ar(i), ar(nx+1-i) )
	814	ENDDO
	815	x_out((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
	816
	817	CALL FFTW_EXECUTE_DFT_C2R( plan_xi, x_out, x_in)
	818	ar(0:nx) = x_in(0:nx)
	819
	820	ENDIF
	821	#endif
	822
[1]	823	ELSEIF ( fft_method == 'system-specific' ) THEN
	824
	825	#if defined( __ibm ) && ! defined( __ibmy_special )
[1106]	826	IF ( forward_fft ) THEN
[1]	827
[1106]	828	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, &
[1]	829	aux2, nau2 )
	830
	831	DO i = 0, (nx+1)/2
	832	ar(i) = work(2*i)
	833	ENDDO
	834	DO i = 1, (nx+1)/2 - 1
	835	ar(nx+1-i) = work(2*i+1)
	836	ENDDO
	837
	838	ELSE
	839
	840	DO i = 0, (nx+1)/2
	841	work(2*i) = ar(i)
	842	ENDDO
	843	DO i = 1, (nx+1)/2 - 1
	844	work(2*i+1) = ar(nx+1-i)
	845	ENDDO
	846	work(1) = 0.0
	847	work(nx+2) = 0.0
	848
[1106]	849	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
[1]	850	aux4, nau2 )
	851
	852	DO i = 0, nx
	853	ar(i) = work(i)
	854	ENDDO
	855
	856	ENDIF
	857	#elif defined( __nec )
[1106]	858	IF ( forward_fft ) THEN
[1]	859
	860	work(0:nx) = ar(0:nx)
	861
[1106]	862	CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
	863
[1]	864	DO i = 0, (nx+1)/2
	865	ar(i) = work(2*i)
	866	ENDDO
	867	DO i = 1, (nx+1)/2 - 1
	868	ar(nx+1-i) = work(2*i+1)
	869	ENDDO
	870
	871	ELSE
	872
	873	DO i = 0, (nx+1)/2
	874	work(2*i) = ar(i)
	875	ENDDO
	876	DO i = 1, (nx+1)/2 - 1
	877	work(2*i+1) = ar(nx+1-i)
	878	ENDDO
	879	work(1) = 0.0
	880	work(nx+2) = 0.0
	881
[1106]	882	CALL ZDFFT( -1, nx+1, sqr_dnx, work, work, trig_xb, work2, 0 )
[1]	883
	884	ar(0:nx) = work(0:nx)
	885
	886	ENDIF
	887	#else
[254]	888	message_string = 'no system-specific fft-call available'
[1106]	889	CALL message( 'fft_x_1d', 'PA0188', 1, 2, 0, 6, 0 )
[1]	890	#endif
	891	ELSE
[274]	892	message_string = 'fft method "' // TRIM( fft_method) // &
	893	'" not available'
[1106]	894	CALL message( 'fft_x_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]	895
	896	ENDIF
	897
[1106]	898	END SUBROUTINE fft_x_1d
[1]	899
[1216]	900	SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
	901	nxr_y_l )
[1]	902
	903	!----------------------------------------------------------------------!
	904	! fft_y !
	905	! !
	906	! Fourier-transformation along y-direction !
[1106]	907	! Version for 2D-decomposition !
[1]	908	! !
	909	! fft_y uses internal algorithms (Singleton or Temperton) or !
	910	! system-specific routines, if they are available !
[1216]	911	! !
	912	! direction: 'forward' or 'backward' !
	913	! ar, ar_tr: 3D data arrays !
	914	! forward: ar: before ar_tr: after transformation !
	915	! backward: ar_tr: before ar: after transfosition !
	916	! !
	917	! In case of non-overlapping transposition/transformation: !
	918	! nxl_y_bound = nxl_y_l = nxl_y !
	919	! nxr_y_bound = nxr_y_l = nxr_y !
	920	! !
	921	! In case of overlapping transposition/transformation !
	922	! - nxl_y_bound and nxr_y_bound have the original values of !
	923	! nxl_y, nxr_y. ar_tr is dimensioned using these values. !
	924	! - nxl_y_l = nxr_y_r. ar is dimensioned with these values, so that !
	925	! transformation is carried out for a 2D-plane only. !
[1]	926	!----------------------------------------------------------------------!
	927
[1106]	928	USE cuda_fft_interfaces
[1153]	929	#if defined( __cuda_fft )
	930	USE ISO_C_BINDING
	931	#endif
[1106]	932
[1]	933	IMPLICIT NONE
	934
	935	CHARACTER (LEN=*) :: direction
[1111]	936	INTEGER :: i, j, jshape(1), k
[1216]	937	INTEGER :: nxl_y_bound, nxl_y_l, nxr_y_bound, nxr_y_l
[1106]	938
	939	LOGICAL :: forward_fft
	940
	941	REAL, DIMENSION(0:ny+2) :: work
	942	REAL, DIMENSION(ny+2) :: work1
	943	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	944	#if defined( __ibm )
	945	REAL, DIMENSION(nau2) :: auy2, auy4
	946	#elif defined( __nec )
	947	REAL, DIMENSION(6*(ny+1)) :: work2
	948	#elif defined( __cuda_fft )
[1111]	949	!$acc declare create( ar_tmp )
[1166]	950	COMPLEX(dpk), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) :: ar_tmp
[1106]	951	#endif
[1216]	952	REAL, DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y) :: ar
	953	REAL, DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) :: ar_tr
[1106]	954
	955	IF ( direction == 'forward' ) THEN
	956	forward_fft = .TRUE.
	957	ELSE
	958	forward_fft = .FALSE.
	959	ENDIF
	960
	961	IF ( fft_method == 'singleton-algorithm' ) THEN
	962
	963	!
	964	!-- Performing the fft with singleton's software works on every system,
	965	!-- since it is part of the model
	966	ALLOCATE( cwork(0:ny) )
	967
	968	IF ( forward_fft ) then
	969
	970	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
	971	!$OMP DO
	972	DO k = nzb_y, nzt_y
[1216]	973	DO i = nxl_y_l, nxr_y_l
[1106]	974
	975	DO j = 0, ny
	976	cwork(j) = CMPLX( ar(j,i,k) )
	977	ENDDO
	978
	979	jshape = SHAPE( cwork )
	980	CALL FFTN( cwork, jshape )
	981
	982	DO j = 0, (ny+1)/2
[1216]	983	ar_tr(j,i,k) = REAL( cwork(j) )
[1106]	984	ENDDO
	985	DO j = 1, (ny+1)/2 - 1
[1216]	986	ar_tr(ny+1-j,i,k) = -AIMAG( cwork(j) )
[1106]	987	ENDDO
	988
	989	ENDDO
	990	ENDDO
	991	!$OMP END PARALLEL
	992
	993	ELSE
	994
	995	!$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
	996	!$OMP DO
	997	DO k = nzb_y, nzt_y
[1216]	998	DO i = nxl_y_l, nxr_y_l
[1106]	999
[1216]	1000	cwork(0) = CMPLX( ar_tr(0,i,k), 0.0 )
[1106]	1001	DO j = 1, (ny+1)/2 - 1
[1216]	1002	cwork(j) = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k) )
	1003	cwork(ny+1-j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
[1106]	1004	ENDDO
[1216]	1005	cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
[1106]	1006
	1007	jshape = SHAPE( cwork )
	1008	CALL FFTN( cwork, jshape, inv = .TRUE. )
	1009
	1010	DO j = 0, ny
	1011	ar(j,i,k) = REAL( cwork(j) )
	1012	ENDDO
	1013
	1014	ENDDO
	1015	ENDDO
	1016	!$OMP END PARALLEL
	1017
	1018	ENDIF
	1019
	1020	DEALLOCATE( cwork )
	1021
	1022	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	1023
	1024	!
	1025	!-- Performing the fft with Temperton's software works on every system,
	1026	!-- since it is part of the model
	1027	IF ( forward_fft ) THEN
	1028
	1029	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1030	!$OMP DO
	1031	DO k = nzb_y, nzt_y
[1216]	1032	DO i = nxl_y_l, nxr_y_l
[1106]	1033
	1034	work(0:ny) = ar(0:ny,i,k)
	1035	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	1036
	1037	DO j = 0, (ny+1)/2
[1216]	1038	ar_tr(j,i,k) = work(2*j)
[1106]	1039	ENDDO
	1040	DO j = 1, (ny+1)/2 - 1
[1216]	1041	ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]	1042	ENDDO
	1043
	1044	ENDDO
	1045	ENDDO
	1046	!$OMP END PARALLEL
	1047
	1048	ELSE
	1049
	1050	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1051	!$OMP DO
	1052	DO k = nzb_y, nzt_y
[1216]	1053	DO i = nxl_y_l, nxr_y_l
[1106]	1054
	1055	DO j = 0, (ny+1)/2
[1216]	1056	work(2*j) = ar_tr(j,i,k)
[1106]	1057	ENDDO
	1058	DO j = 1, (ny+1)/2 - 1
[1216]	1059	work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]	1060	ENDDO
	1061	work(1) = 0.0
	1062	work(ny+2) = 0.0
	1063
	1064	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	1065	ar(0:ny,i,k) = work(0:ny)
	1066
	1067	ENDDO
	1068	ENDDO
	1069	!$OMP END PARALLEL
	1070
	1071	ENDIF
	1072
[1210]	1073	ELSEIF ( fft_method == 'fftw' ) THEN
	1074
	1075	#if defined( __fftw )
	1076	IF ( forward_fft ) THEN
	1077
	1078	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1079	!$OMP DO
	1080	DO k = nzb_y, nzt_y
[1216]	1081	DO i = nxl_y_l, nxr_y_l
[1210]	1082
	1083	y_in(0:ny) = ar(0:ny,i,k)
	1084	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
	1085
	1086	DO j = 0, (ny+1)/2
[1216]	1087	ar_tr(j,i,k) = REAL( y_out(j) ) / (ny+1)
[1210]	1088	ENDDO
	1089	DO j = 1, (ny+1)/2 - 1
[1216]	1090	ar_tr(ny+1-j,i,k) = AIMAG( y_out(j) ) / (ny+1)
[1210]	1091	ENDDO
	1092
	1093	ENDDO
	1094	ENDDO
	1095	!$OMP END PARALLEL
	1096
	1097	ELSE
	1098
	1099	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1100	!$OMP DO
	1101	DO k = nzb_y, nzt_y
[1216]	1102	DO i = nxl_y_l, nxr_y_l
[1210]	1103
[1216]	1104	y_out(0) = CMPLX( ar_tr(0,i,k), 0.0 )
[1210]	1105	DO j = 1, (ny+1)/2 - 1
[1216]	1106	y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k) )
[1210]	1107	ENDDO
[1216]	1108	y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0 )
[1210]	1109
	1110	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
	1111	ar(0:ny,i,k) = y_in(0:ny)
	1112
	1113	ENDDO
	1114	ENDDO
	1115	!$OMP END PARALLEL
	1116
	1117	ENDIF
	1118	#endif
	1119
[1106]	1120	ELSEIF ( fft_method == 'system-specific' ) THEN
	1121
	1122	#if defined( __ibm ) && ! defined( __ibmy_special )
	1123	IF ( forward_fft) THEN
	1124
	1125	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1126	!$OMP DO
	1127	DO k = nzb_y, nzt_y
[1216]	1128	DO i = nxl_y_l, nxr_y_l
[1106]	1129
	1130	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
	1131	auy2, nau2 )
	1132
	1133	DO j = 0, (ny+1)/2
[1216]	1134	ar_tr(j,i,k) = work(2*j)
[1106]	1135	ENDDO
	1136	DO j = 1, (ny+1)/2 - 1
[1216]	1137	ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]	1138	ENDDO
	1139
	1140	ENDDO
	1141	ENDDO
	1142	!$OMP END PARALLEL
	1143
	1144	ELSE
	1145
	1146	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1147	!$OMP DO
	1148	DO k = nzb_y, nzt_y
[1216]	1149	DO i = nxl_y_l, nxr_y_l
[1106]	1150
	1151	DO j = 0, (ny+1)/2
[1216]	1152	work(2*j) = ar_tr(j,i,k)
[1106]	1153	ENDDO
	1154	DO j = 1, (ny+1)/2 - 1
[1216]	1155	work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]	1156	ENDDO
	1157	work(1) = 0.0
	1158	work(ny+2) = 0.0
	1159
	1160	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
	1161	auy4, nau2 )
	1162
	1163	DO j = 0, ny
	1164	ar(j,i,k) = work(j)
	1165	ENDDO
	1166
	1167	ENDDO
	1168	ENDDO
	1169	!$OMP END PARALLEL
	1170
	1171	ENDIF
	1172	#elif defined( __nec )
	1173	IF ( forward_fft ) THEN
	1174
	1175	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1176	!$OMP DO
	1177	DO k = nzb_y, nzt_y
[1216]	1178	DO i = nxl_y_l, nxr_y_l
[1106]	1179
	1180	work(0:ny) = ar(0:ny,i,k)
	1181
	1182	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
	1183
	1184	DO j = 0, (ny+1)/2
[1216]	1185	ar_tr(j,i,k) = work(2*j)
[1106]	1186	ENDDO
	1187	DO j = 1, (ny+1)/2 - 1
[1216]	1188	ar_tr(ny+1-j,i,k) = work(2*j+1)
[1106]	1189	ENDDO
	1190
	1191	ENDDO
	1192	ENDDO
	1193	!$END OMP PARALLEL
	1194
	1195	ELSE
	1196
	1197	!$OMP PARALLEL PRIVATE ( work, i, j, k )
	1198	!$OMP DO
	1199	DO k = nzb_y, nzt_y
[1216]	1200	DO i = nxl_y_l, nxr_y_l
[1106]	1201
	1202	DO j = 0, (ny+1)/2
[1216]	1203	work(2*j) = ar_tr(j,i,k)
[1106]	1204	ENDDO
	1205	DO j = 1, (ny+1)/2 - 1
[1216]	1206	work(2*j+1) = ar_tr(ny+1-j,i,k)
[1106]	1207	ENDDO
	1208	work(1) = 0.0
	1209	work(ny+2) = 0.0
	1210
	1211	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
	1212
	1213	ar(0:ny,i,k) = work(0:ny)
	1214
	1215	ENDDO
	1216	ENDDO
	1217	!$OMP END PARALLEL
	1218
	1219	ENDIF
	1220	#elif defined( __cuda_fft )
	1221
	1222	IF ( forward_fft ) THEN
	1223
[1111]	1224	!$acc data present( ar )
	1225	CALL CUFFTEXECD2Z( plan_yf, ar, ar_tmp )
[1106]	1226
[1111]	1227	!$acc kernels
	1228	!$acc loop
[1106]	1229	DO k = nzb_y, nzt_y
	1230	DO i = nxl_y, nxr_y
	1231
[1111]	1232	!$acc loop vector( 32 )
[1106]	1233	DO j = 0, (ny+1)/2
[1111]	1234	ar(j,i,k) = REAL( ar_tmp(j,i,k) ) * dny
[1106]	1235	ENDDO
	1236
[1111]	1237	!$acc loop vector( 32 )
[1106]	1238	DO j = 1, (ny+1)/2 - 1
[1111]	1239	ar(ny+1-j,i,k) = AIMAG( ar_tmp(j,i,k) ) * dny
[1106]	1240	ENDDO
	1241
	1242	ENDDO
	1243	ENDDO
[1111]	1244	!$acc end kernels
	1245	!$acc end data
[1106]	1246
	1247	ELSE
	1248
[1111]	1249	!$acc data present( ar )
	1250	!$acc kernels
	1251	!$acc loop
[1106]	1252	DO k = nzb_y, nzt_y
	1253	DO i = nxl_y, nxr_y
	1254
[1111]	1255	ar_tmp(0,i,k) = CMPLX( ar(0,i,k), 0.0 )
[1106]	1256
[1111]	1257	!$acc loop vector( 32 )
[1106]	1258	DO j = 1, (ny+1)/2 - 1
[1111]	1259	ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k) )
[1106]	1260	ENDDO
[1111]	1261	ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0 )
[1106]	1262
	1263	ENDDO
	1264	ENDDO
[1111]	1265	!$acc end kernels
[1106]	1266
[1111]	1267	CALL CUFFTEXECZ2D( plan_yi, ar_tmp, ar )
	1268	!$acc end data
[1106]	1269
	1270	ENDIF
	1271
	1272	#else
	1273	message_string = 'no system-specific fft-call available'
	1274	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
	1275	#endif
	1276
	1277	ELSE
	1278
	1279	message_string = 'fft method "' // TRIM( fft_method) // &
	1280	'" not available'
	1281	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
	1282
	1283	ENDIF
	1284
	1285	END SUBROUTINE fft_y
	1286
	1287	SUBROUTINE fft_y_1d( ar, direction )
	1288
	1289	!----------------------------------------------------------------------!
	1290	! fft_y_1d !
	1291	! !
	1292	! Fourier-transformation along y-direction !
	1293	! Version for 1D-decomposition !
	1294	! !
	1295	! fft_y uses internal algorithms (Singleton or Temperton) or !
	1296	! system-specific routines, if they are available !
	1297	!----------------------------------------------------------------------!
	1298
	1299	IMPLICIT NONE
	1300
	1301	CHARACTER (LEN=*) :: direction
[1]	1302	INTEGER :: j, jshape(1)
	1303
[1106]	1304	LOGICAL :: forward_fft
	1305
[1]	1306	REAL, DIMENSION(0:ny) :: ar
	1307	REAL, DIMENSION(0:ny+2) :: work
	1308	REAL, DIMENSION(ny+2) :: work1
	1309	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	1310	#if defined( __ibm )
	1311	REAL, DIMENSION(nau2) :: auy2, auy4
	1312	#elif defined( __nec )
	1313	REAL, DIMENSION(6*(ny+1)) :: work2
	1314	#endif
	1315
[1106]	1316	IF ( direction == 'forward' ) THEN
	1317	forward_fft = .TRUE.
	1318	ELSE
	1319	forward_fft = .FALSE.
	1320	ENDIF
	1321
[1]	1322	IF ( fft_method == 'singleton-algorithm' ) THEN
	1323
	1324	!
	1325	!-- Performing the fft with singleton's software works on every system,
	1326	!-- since it is part of the model
	1327	ALLOCATE( cwork(0:ny) )
	1328
[1106]	1329	IF ( forward_fft ) THEN
[1]	1330
	1331	DO j = 0, ny
	1332	cwork(j) = CMPLX( ar(j) )
	1333	ENDDO
	1334
	1335	jshape = SHAPE( cwork )
	1336	CALL FFTN( cwork, jshape )
	1337
	1338	DO j = 0, (ny+1)/2
	1339	ar(j) = REAL( cwork(j) )
	1340	ENDDO
	1341	DO j = 1, (ny+1)/2 - 1
	1342	ar(ny+1-j) = -AIMAG( cwork(j) )
	1343	ENDDO
	1344
	1345	ELSE
	1346
	1347	cwork(0) = CMPLX( ar(0), 0.0 )
	1348	DO j = 1, (ny+1)/2 - 1
	1349	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
	1350	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
	1351	ENDDO
	1352	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
	1353
	1354	jshape = SHAPE( cwork )
	1355	CALL FFTN( cwork, jshape, inv = .TRUE. )
	1356
	1357	DO j = 0, ny
	1358	ar(j) = REAL( cwork(j) )
	1359	ENDDO
	1360
	1361	ENDIF
	1362
	1363	DEALLOCATE( cwork )
	1364
	1365	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	1366
	1367	!
	1368	!-- Performing the fft with Temperton's software works on every system,
	1369	!-- since it is part of the model
[1106]	1370	IF ( forward_fft ) THEN
[1]	1371
	1372	work(0:ny) = ar
	1373	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	1374
	1375	DO j = 0, (ny+1)/2
	1376	ar(j) = work(2*j)
	1377	ENDDO
	1378	DO j = 1, (ny+1)/2 - 1
	1379	ar(ny+1-j) = work(2*j+1)
	1380	ENDDO
	1381
	1382	ELSE
	1383
	1384	DO j = 0, (ny+1)/2
	1385	work(2*j) = ar(j)
	1386	ENDDO
	1387	DO j = 1, (ny+1)/2 - 1
	1388	work(2*j+1) = ar(ny+1-j)
	1389	ENDDO
	1390	work(1) = 0.0
	1391	work(ny+2) = 0.0
	1392
	1393	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	1394	ar = work(0:ny)
	1395
	1396	ENDIF
	1397
[1216]	1398	ELSEIF ( fft_method == 'fftw' ) THEN
	1399
	1400	#if defined( __fftw )
	1401	IF ( forward_fft ) THEN
	1402
	1403	y_in(0:ny) = ar(0:ny)
	1404	CALL FFTW_EXECUTE_DFT_R2C( plan_yf, y_in, y_out )
	1405
	1406	DO j = 0, (ny+1)/2
	1407	ar(j) = REAL( y_out(j) ) / (ny+1)
	1408	ENDDO
	1409	DO j = 1, (ny+1)/2 - 1
	1410	ar(ny+1-j) = AIMAG( y_out(j) ) / (ny+1)
	1411	ENDDO
	1412
	1413	ELSE
	1414
	1415	y_out(0) = CMPLX( ar(0), 0.0 )
	1416	DO j = 1, (ny+1)/2 - 1
	1417	y_out(j) = CMPLX( ar(j), ar(ny+1-j) )
	1418	ENDDO
	1419	y_out((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
	1420
	1421	CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
	1422	ar(0:ny) = y_in(0:ny)
	1423
	1424	ENDIF
	1425	#endif
	1426
[1]	1427	ELSEIF ( fft_method == 'system-specific' ) THEN
	1428
	1429	#if defined( __ibm ) && ! defined( __ibmy_special )
[1106]	1430	IF ( forward_fft ) THEN
[1]	1431
[1106]	1432	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, &
[1]	1433	auy2, nau2 )
	1434
	1435	DO j = 0, (ny+1)/2
	1436	ar(j) = work(2*j)
	1437	ENDDO
	1438	DO j = 1, (ny+1)/2 - 1
	1439	ar(ny+1-j) = work(2*j+1)
	1440	ENDDO
	1441
	1442	ELSE
	1443
	1444	DO j = 0, (ny+1)/2
	1445	work(2*j) = ar(j)
	1446	ENDDO
	1447	DO j = 1, (ny+1)/2 - 1
	1448	work(2*j+1) = ar(ny+1-j)
	1449	ENDDO
	1450	work(1) = 0.0
	1451	work(ny+2) = 0.0
	1452
[1106]	1453	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
[1]	1454	auy4, nau2 )
	1455
	1456	DO j = 0, ny
	1457	ar(j) = work(j)
	1458	ENDDO
	1459
	1460	ENDIF
	1461	#elif defined( __nec )
[1106]	1462	IF ( forward_fft ) THEN
[1]	1463
	1464	work(0:ny) = ar(0:ny)
	1465
[1106]	1466	CALL DZFFT( 1, ny+1, sqr_dny, work, work, trig_yf, work2, 0 )
[1]	1467
	1468	DO j = 0, (ny+1)/2
	1469	ar(j) = work(2*j)
	1470	ENDDO
	1471	DO j = 1, (ny+1)/2 - 1
	1472	ar(ny+1-j) = work(2*j+1)
	1473	ENDDO
	1474
	1475	ELSE
	1476
	1477	DO j = 0, (ny+1)/2
	1478	work(2*j) = ar(j)
	1479	ENDDO
	1480	DO j = 1, (ny+1)/2 - 1
	1481	work(2*j+1) = ar(ny+1-j)
	1482	ENDDO
	1483	work(1) = 0.0
	1484	work(ny+2) = 0.0
	1485
[1106]	1486	CALL ZDFFT( -1, ny+1, sqr_dny, work, work, trig_yb, work2, 0 )
[1]	1487
	1488	ar(0:ny) = work(0:ny)
	1489
	1490	ENDIF
	1491	#else
[254]	1492	message_string = 'no system-specific fft-call available'
[1106]	1493	CALL message( 'fft_y_1d', 'PA0188', 1, 2, 0, 6, 0 )
[254]	1494
[1]	1495	#endif
	1496
	1497	ELSE
	1498
[274]	1499	message_string = 'fft method "' // TRIM( fft_method) // &
	1500	'" not available'
[1106]	1501	CALL message( 'fft_y_1d', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1502
	1503	ENDIF
	1504
[1106]	1505	END SUBROUTINE fft_y_1d
[1]	1506
	1507	SUBROUTINE fft_x_m( ar, direction )
	1508
	1509	!----------------------------------------------------------------------!
	1510	! fft_x_m !
	1511	! !
	1512	! Fourier-transformation along x-direction !
	1513	! Version for 1d domain decomposition !
	1514	! using multiple 1D FFT from Math Keisan on NEC !
	1515	! or Temperton-algorithm !
	1516	! (no singleton-algorithm on NEC because it does not vectorize) !
	1517	! !
	1518	!----------------------------------------------------------------------!
	1519
	1520	IMPLICIT NONE
	1521
	1522	CHARACTER (LEN=*) :: direction
[1092]	1523	INTEGER :: i, k, siza
[1]	1524
	1525	REAL, DIMENSION(0:nx,nz) :: ar
	1526	REAL, DIMENSION(0:nx+3,nz+1) :: ai
	1527	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
	1528	#if defined( __nec )
[1092]	1529	INTEGER :: sizw
[1]	1530	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
	1531	#endif
	1532
	1533	IF ( fft_method == 'temperton-algorithm' ) THEN
	1534
	1535	siza = SIZE( ai, 1 )
	1536
	1537	IF ( direction == 'forward') THEN
	1538
	1539	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	1540	ai(nx+1:,:) = 0.0
	1541
	1542	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
	1543
	1544	DO k = 1, nz
	1545	DO i = 0, (nx+1)/2
	1546	ar(i,k) = ai(2*i,k)
	1547	ENDDO
	1548	DO i = 1, (nx+1)/2 - 1
	1549	ar(nx+1-i,k) = ai(2*i+1,k)
	1550	ENDDO
	1551	ENDDO
	1552
	1553	ELSE
	1554
	1555	DO k = 1, nz
	1556	DO i = 0, (nx+1)/2
	1557	ai(2*i,k) = ar(i,k)
	1558	ENDDO
	1559	DO i = 1, (nx+1)/2 - 1
	1560	ai(2*i+1,k) = ar(nx+1-i,k)
	1561	ENDDO
	1562	ai(1,k) = 0.0
	1563	ai(nx+2,k) = 0.0
	1564	ENDDO
	1565
	1566	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
	1567
	1568	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	1569
	1570	ENDIF
	1571
	1572	ELSEIF ( fft_method == 'system-specific' ) THEN
	1573
	1574	#if defined( __nec )
	1575	siza = SIZE( ai, 1 )
	1576	sizw = SIZE( work, 1 )
	1577
	1578	IF ( direction == 'forward') THEN
	1579
	1580	!
	1581	!-- Tables are initialized once more. This call should not be
	1582	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1583	CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]	1584	trig_xf, work1, 0 )
	1585
	1586	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	1587	IF ( nz1 > nz ) THEN
	1588	ai(:,nz1) = 0.0
	1589	ENDIF
	1590
[1106]	1591	CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, &
[1]	1592	trig_xf, work1, 0 )
	1593
	1594	DO k = 1, nz
	1595	DO i = 0, (nx+1)/2
	1596	ar(i,k) = REAL( work(i+1,k) )
	1597	ENDDO
	1598	DO i = 1, (nx+1)/2 - 1
	1599	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
	1600	ENDDO
	1601	ENDDO
	1602
	1603	ELSE
	1604
	1605	!
	1606	!-- Tables are initialized once more. This call should not be
	1607	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1608	CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, &
[1]	1609	trig_xb, work1, 0 )
	1610
	1611	IF ( nz1 > nz ) THEN
	1612	work(:,nz1) = 0.0
	1613	ENDIF
	1614	DO k = 1, nz
	1615	work(1,k) = CMPLX( ar(0,k), 0.0 )
	1616	DO i = 1, (nx+1)/2 - 1
	1617	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
	1618	ENDDO
	1619	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
	1620	ENDDO
	1621
[1106]	1622	CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
[1]	1623	trig_xb, work1, 0 )
	1624
	1625	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	1626
	1627	ENDIF
	1628
	1629	#else
[254]	1630	message_string = 'no system-specific fft-call available'
	1631	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	1632	#endif
	1633
	1634	ELSE
	1635
[274]	1636	message_string = 'fft method "' // TRIM( fft_method) // &
	1637	'" not available'
[254]	1638	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1639
	1640	ENDIF
	1641
	1642	END SUBROUTINE fft_x_m
	1643
	1644	SUBROUTINE fft_y_m( ar, ny1, direction )
	1645
	1646	!----------------------------------------------------------------------!
	1647	! fft_y_m !
	1648	! !
	1649	! Fourier-transformation along y-direction !
	1650	! Version for 1d domain decomposition !
	1651	! using multiple 1D FFT from Math Keisan on NEC !
	1652	! or Temperton-algorithm !
	1653	! (no singleton-algorithm on NEC because it does not vectorize) !
	1654	! !
	1655	!----------------------------------------------------------------------!
	1656
	1657	IMPLICIT NONE
	1658
	1659	CHARACTER (LEN=*) :: direction
[1092]	1660	INTEGER :: j, k, ny1, siza
[1]	1661
	1662	REAL, DIMENSION(0:ny1,nz) :: ar
	1663	REAL, DIMENSION(0:ny+3,nz+1) :: ai
	1664	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
	1665	#if defined( __nec )
[1092]	1666	INTEGER :: sizw
[1]	1667	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
	1668	#endif
	1669
	1670	IF ( fft_method == 'temperton-algorithm' ) THEN
	1671
	1672	siza = SIZE( ai, 1 )
	1673
	1674	IF ( direction == 'forward') THEN
	1675
	1676	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	1677	ai(ny+1:,:) = 0.0
	1678
	1679	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
	1680
	1681	DO k = 1, nz
	1682	DO j = 0, (ny+1)/2
	1683	ar(j,k) = ai(2*j,k)
	1684	ENDDO
	1685	DO j = 1, (ny+1)/2 - 1
	1686	ar(ny+1-j,k) = ai(2*j+1,k)
	1687	ENDDO
	1688	ENDDO
	1689
	1690	ELSE
	1691
	1692	DO k = 1, nz
	1693	DO j = 0, (ny+1)/2
	1694	ai(2*j,k) = ar(j,k)
	1695	ENDDO
	1696	DO j = 1, (ny+1)/2 - 1
	1697	ai(2*j+1,k) = ar(ny+1-j,k)
	1698	ENDDO
	1699	ai(1,k) = 0.0
	1700	ai(ny+2,k) = 0.0
	1701	ENDDO
	1702
	1703	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
	1704
	1705	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	1706
	1707	ENDIF
	1708
	1709	ELSEIF ( fft_method == 'system-specific' ) THEN
	1710
	1711	#if defined( __nec )
	1712	siza = SIZE( ai, 1 )
	1713	sizw = SIZE( work, 1 )
	1714
	1715	IF ( direction == 'forward') THEN
	1716
	1717	!
	1718	!-- Tables are initialized once more. This call should not be
	1719	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1720	CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]	1721	trig_yf, work1, 0 )
	1722
	1723	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	1724	IF ( nz1 > nz ) THEN
	1725	ai(:,nz1) = 0.0
	1726	ENDIF
	1727
[1106]	1728	CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
[1]	1729	trig_yf, work1, 0 )
	1730
	1731	DO k = 1, nz
	1732	DO j = 0, (ny+1)/2
	1733	ar(j,k) = REAL( work(j+1,k) )
	1734	ENDDO
	1735	DO j = 1, (ny+1)/2 - 1
	1736	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
	1737	ENDDO
	1738	ENDDO
	1739
	1740	ELSE
	1741
	1742	!
	1743	!-- Tables are initialized once more. This call should not be
	1744	!-- necessary, but otherwise program aborts in asymmetric case
[1106]	1745	CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
[1]	1746	trig_yb, work1, 0 )
	1747
	1748	IF ( nz1 > nz ) THEN
	1749	work(:,nz1) = 0.0
	1750	ENDIF
	1751	DO k = 1, nz
	1752	work(1,k) = CMPLX( ar(0,k), 0.0 )
	1753	DO j = 1, (ny+1)/2 - 1
	1754	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
	1755	ENDDO
	1756	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
	1757	ENDDO
	1758
[1106]	1759	CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
[1]	1760	trig_yb, work1, 0 )
	1761
	1762	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	1763
	1764	ENDIF
	1765
	1766	#else
[254]	1767	message_string = 'no system-specific fft-call available'
	1768	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	1769	#endif
	1770
	1771	ELSE
[254]	1772
[274]	1773	message_string = 'fft method "' // TRIM( fft_method) // &
	1774	'" not available'
[254]	1775	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	1776
	1777	ENDIF
	1778
	1779	END SUBROUTINE fft_y_m
	1780
[1106]	1781
[1]	1782	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |