Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |