Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |