Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |