Home

Context Navigation

fft_xy.f90 @ 1111

Last change on this file since 1111 was 1111, checked in by raasch, 11 years ago

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:

resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:

bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: