Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |