Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |