Home

Context Navigation

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

Last change on this file since 1037 was 1037, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 25.5 KB

Rev	Line
[1]	1	MODULE fft_xy
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[392]	22	!
[1]	23	!
	24	! Former revisions:
	25	! -----------------
[3]	26	! $Id: fft_xy.f90 1037 2012-10-22 14:10:22Z raasch $
[392]	27	!
[1037]	28	! 1036 2012-10-22 13:43:42Z raasch
	29	! code put under GPL (PALM 3.9)
	30	!
[392]	31	! 274 2009-03-26 15:11:21Z heinze
	32	! Output of messages replaced by message handling routine.
	33	!
	34	! Feb. 2007
[3]	35	! RCS Log replace by Id keyword, revision history cleaned up
	36	!
[1]	37	! Revision 1.4 2006/03/28 12:27:09 raasch
	38	! Stop when system-specific fft is selected on NEC. For unknown reasons this
	39	! causes a program abort during first allocation in init_grid.
	40	!
	41	! Revision 1.2 2004/04/30 11:44:27 raasch
	42	! Module renamed from fft_for_1d_decomp to fft_xy, 1d-routines renamed to
	43	! fft_x and fft_y,
	44	! function FFT replaced by subroutine FFTN due to problems with 64-bit
	45	! mode on ibm,
	46	! shape of array cwork is explicitly stored in ishape/jshape and handled
	47	! to routine FFTN instead of shape-function (due to compiler error on
	48	! decalpha),
	49	! non vectorized FFT for nec included
	50	!
	51	! Revision 1.1 2002/06/11 13:00:49 raasch
	52	! Initial revision
	53	!
	54	!
	55	! Description:
	56	! ------------
	57	! Fast Fourier transformation along x and y for 1d domain decomposition along x.
	58	! Original version: Klaus Ketelsen (May 2002)
	59	!------------------------------------------------------------------------------!
	60
	61	USE array_kind
	62	USE control_parameters
	63	USE indices
	64	USE singleton
	65	USE temperton_fft
	66
	67	IMPLICIT NONE
	68
	69	PRIVATE
	70	PUBLIC fft_x, fft_y, fft_init, fft_x_m, fft_y_m
	71
	72	INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ifax_x, ifax_y
	73
	74	LOGICAL, SAVE :: init_fft = .FALSE.
	75
	76	REAL, SAVE :: sqr_nx, sqr_ny
	77	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trigs_x, trigs_y
	78
	79	#if defined( __ibm )
	80	INTEGER, PARAMETER :: nau1 = 20000, nau2 = 22000
	81	!
	82	!-- The following working arrays contain tables and have to be "save" and
	83	!-- shared in OpenMP sense
	84	REAL, DIMENSION(nau1), SAVE :: aux1, auy1, aux3, auy3
	85	#elif defined( __nec )
	86	INTEGER, SAVE :: nz1
	87	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: trig_xb, trig_xf, trig_yb, &
	88	trig_yf
	89	#endif
	90
	91	!
	92	!-- Public interfaces
	93	INTERFACE fft_init
	94	MODULE PROCEDURE fft_init
	95	END INTERFACE fft_init
	96
	97	INTERFACE fft_x
	98	MODULE PROCEDURE fft_x
	99	END INTERFACE fft_x
	100
	101	INTERFACE fft_y
	102	MODULE PROCEDURE fft_y
	103	END INTERFACE fft_y
	104
	105	INTERFACE fft_x_m
	106	MODULE PROCEDURE fft_x_m
	107	END INTERFACE fft_x_m
	108
	109	INTERFACE fft_y_m
	110	MODULE PROCEDURE fft_y_m
	111	END INTERFACE fft_y_m
	112
	113	CONTAINS
	114
	115
	116	SUBROUTINE fft_init
	117
	118	IMPLICIT NONE
	119
	120	!
	121	!-- The following temporary working arrays have to be on stack or private
	122	!-- in OpenMP sense
	123	#if defined( __ibm )
	124	REAL, DIMENSION(0:nx+2) :: workx
	125	REAL, DIMENSION(0:ny+2) :: worky
	126	REAL, DIMENSION(nau2) :: aux2, auy2, aux4, auy4
	127	#elif defined( __nec )
	128	REAL, DIMENSION(0:nx+3,nz+1) :: work_x
	129	REAL, DIMENSION(0:ny+3,nz+1) :: work_y
	130	REAL, DIMENSION(6*(nx+3),nz+1) :: workx
	131	REAL, DIMENSION(6*(ny+3),nz+1) :: worky
	132	#endif
	133
	134	!
	135	!-- Return, if already called
	136	IF ( init_fft ) THEN
	137	RETURN
	138	ELSE
	139	init_fft = .TRUE.
	140	ENDIF
	141
	142	IF ( fft_method == 'system-specific' ) THEN
	143
	144	sqr_nx = SQRT( 1.0 / ( nx + 1.0 ) )
	145	sqr_ny = SQRT( 1.0 / ( ny + 1.0 ) )
	146	#if defined( __ibm ) && ! defined( __ibmy_special )
	147	!
	148	!-- Initialize tables for fft along x
	149	CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
	150	aux2, nau2 )
	151	CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
	152	aux4, nau2 )
	153	!
	154	!-- Initialize tables for fft along y
	155	CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
	156	auy2, nau2 )
	157	CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
	158	auy4, nau2 )
	159	#elif defined( __nec )
[254]	160	message_string = 'fft method "' // TRIM( fft_method) // &
	161	'" currently does not work on NEC'
	162	CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
[1]	163
	164	ALLOCATE( trig_xb(2(nx+1)), trig_xf(2(nx+1)), &
	165	trig_yb(2(ny+1)), trig_yf(2(ny+1)) )
	166
	167	work_x = 0.0
	168	work_y = 0.0
	169	nz1 = nz + MOD( nz+1, 2 ) ! odd nz slows down fft significantly
	170	! when using the NEC ffts
	171
	172	!
	173	!-- Initialize tables for fft along x (non-vector and vector case (M))
	174	CALL DZFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xf, workx, 0 )
	175	CALL ZDFFT( 0, nx+1, sqr_nx, work_x, work_x, trig_xb, workx, 0 )
	176	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
	177	trig_xf, workx, 0 )
	178	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work_x, nx+4, work_x, nx+4, &
	179	trig_xb, workx, 0 )
	180	!
	181	!-- Initialize tables for fft along y (non-vector and vector case (M))
	182	CALL DZFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yf, worky, 0 )
	183	CALL ZDFFT( 0, ny+1, sqr_ny, work_y, work_y, trig_yb, worky, 0 )
	184	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
	185	trig_yf, worky, 0 )
	186	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work_y, ny+4, work_y, ny+4, &
	187	trig_yb, worky, 0 )
	188	#else
[254]	189	message_string = 'no system-specific fft-call available'
	190	CALL message( 'fft_init', 'PA0188', 1, 2, 0, 6, 0 )
[1]	191	#endif
	192	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	193	!
	194	!-- Temperton-algorithm
	195	!-- Initialize tables for fft along x and y
	196	ALLOCATE( ifax_x(nx+1), ifax_y(ny+1), trigs_x(nx+1), trigs_y(ny+1) )
	197
	198	CALL set99( trigs_x, ifax_x, nx+1 )
	199	CALL set99( trigs_y, ifax_y, ny+1 )
	200
	201	ELSEIF ( fft_method == 'singleton-algorithm' ) THEN
	202
	203	CONTINUE
	204
	205	ELSE
	206
[254]	207	message_string = 'fft method "' // TRIM( fft_method) // &
	208	'" not available'
	209	CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
[1]	210	ENDIF
	211
	212	END SUBROUTINE fft_init
	213
	214
	215	SUBROUTINE fft_x( ar, direction )
	216
	217	!----------------------------------------------------------------------!
	218	! fft_x !
	219	! !
	220	! Fourier-transformation along x-direction !
	221	! !
	222	! fft_x uses internal algorithms (Singleton or Temperton) or !
	223	! system-specific routines, if they are available !
	224	!----------------------------------------------------------------------!
	225
	226	IMPLICIT NONE
	227
	228	CHARACTER (LEN=*) :: direction
	229	INTEGER :: i, ishape(1)
	230
	231	!kk REAL, DIMENSION(:) :: ar !kk Does NOT work (Bug??)
	232	REAL, DIMENSION(0:nx) :: ar
	233	REAL, DIMENSION(0:nx+2) :: work
	234	REAL, DIMENSION(nx+2) :: work1
	235	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	236	#if defined( __ibm )
	237	REAL, DIMENSION(nau2) :: aux2, aux4
	238	#elif defined( __nec )
	239	REAL, DIMENSION(6*(nx+1)) :: work2
	240	#endif
	241
	242	IF ( fft_method == 'singleton-algorithm' ) THEN
	243
	244	!
	245	!-- Performing the fft with singleton's software works on every system,
	246	!-- since it is part of the model
	247	ALLOCATE( cwork(0:nx) )
	248
	249	IF ( direction == 'forward') then
	250
	251	DO i = 0, nx
	252	cwork(i) = CMPLX( ar(i) )
	253	ENDDO
	254	ishape = SHAPE( cwork )
	255	CALL FFTN( cwork, ishape )
	256
	257	DO i = 0, (nx+1)/2
	258	ar(i) = REAL( cwork(i) )
	259	ENDDO
	260	DO i = 1, (nx+1)/2 - 1
	261	ar(nx+1-i) = -AIMAG( cwork(i) )
	262	ENDDO
	263
	264	ELSE
	265
	266	cwork(0) = CMPLX( ar(0), 0.0 )
	267	DO i = 1, (nx+1)/2 - 1
	268	cwork(i) = CMPLX( ar(i), -ar(nx+1-i) )
	269	cwork(nx+1-i) = CMPLX( ar(i), ar(nx+1-i) )
	270	ENDDO
	271	cwork((nx+1)/2) = CMPLX( ar((nx+1)/2), 0.0 )
	272
	273	ishape = SHAPE( cwork )
	274	CALL FFTN( cwork, ishape, inv = .TRUE. )
	275
	276	DO i = 0, nx
	277	ar(i) = REAL( cwork(i) )
	278	ENDDO
	279
	280	ENDIF
	281
	282	DEALLOCATE( cwork )
	283
	284	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	285
	286	!
	287	!-- Performing the fft with Temperton's software works on every system,
	288	!-- since it is part of the model
	289	IF ( direction == 'forward' ) THEN
	290
	291	work(0:nx) = ar
	292	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
	293
	294	DO i = 0, (nx+1)/2
	295	ar(i) = work(2*i)
	296	ENDDO
	297	DO i = 1, (nx+1)/2 - 1
	298	ar(nx+1-i) = work(2*i+1)
	299	ENDDO
	300
	301	ELSE
	302
	303	DO i = 0, (nx+1)/2
	304	work(2*i) = ar(i)
	305	ENDDO
	306	DO i = 1, (nx+1)/2 - 1
	307	work(2*i+1) = ar(nx+1-i)
	308	ENDDO
	309	work(1) = 0.0
	310	work(nx+2) = 0.0
	311
	312	CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
	313	ar = work(0:nx)
	314
	315	ENDIF
	316
	317	ELSEIF ( fft_method == 'system-specific' ) THEN
	318
	319	#if defined( __ibm ) && ! defined( __ibmy_special )
	320	IF ( direction == 'forward' ) THEN
	321
	322	CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_nx, aux1, nau1, &
	323	aux2, nau2 )
	324
	325	DO i = 0, (nx+1)/2
	326	ar(i) = work(2*i)
	327	ENDDO
	328	DO i = 1, (nx+1)/2 - 1
	329	ar(nx+1-i) = work(2*i+1)
	330	ENDDO
	331
	332	ELSE
	333
	334	DO i = 0, (nx+1)/2
	335	work(2*i) = ar(i)
	336	ENDDO
	337	DO i = 1, (nx+1)/2 - 1
	338	work(2*i+1) = ar(nx+1-i)
	339	ENDDO
	340	work(1) = 0.0
	341	work(nx+2) = 0.0
	342
	343	CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_nx, aux3, nau1, &
	344	aux4, nau2 )
	345
	346	DO i = 0, nx
	347	ar(i) = work(i)
	348	ENDDO
	349
	350	ENDIF
	351	#elif defined( __nec )
	352	IF ( direction == 'forward' ) THEN
	353
	354	work(0:nx) = ar(0:nx)
	355
	356	CALL DZFFT( 1, nx+1, sqr_nx, work, work, trig_xf, work2, 0 )
	357
	358	DO i = 0, (nx+1)/2
	359	ar(i) = work(2*i)
	360	ENDDO
	361	DO i = 1, (nx+1)/2 - 1
	362	ar(nx+1-i) = work(2*i+1)
	363	ENDDO
	364
	365	ELSE
	366
	367	DO i = 0, (nx+1)/2
	368	work(2*i) = ar(i)
	369	ENDDO
	370	DO i = 1, (nx+1)/2 - 1
	371	work(2*i+1) = ar(nx+1-i)
	372	ENDDO
	373	work(1) = 0.0
	374	work(nx+2) = 0.0
	375
	376	CALL ZDFFT( -1, nx+1, sqr_nx, work, work, trig_xb, work2, 0 )
	377
	378	ar(0:nx) = work(0:nx)
	379
	380	ENDIF
	381	#else
[254]	382	message_string = 'no system-specific fft-call available'
	383	CALL message( 'fft_x', 'PA0188', 1, 2, 0, 6, 0 )
[1]	384	#endif
	385	ELSE
[274]	386	message_string = 'fft method "' // TRIM( fft_method) // &
	387	'" not available'
[254]	388	CALL message( 'fft_x', 'PA0189', 1, 2, 0, 6, 0 )
[1]	389
	390	ENDIF
	391
	392	END SUBROUTINE fft_x
	393
	394	SUBROUTINE fft_y( ar, direction )
	395
	396	!----------------------------------------------------------------------!
	397	! fft_y !
	398	! !
	399	! Fourier-transformation along y-direction !
	400	! !
	401	! fft_y uses internal algorithms (Singleton or Temperton) or !
	402	! system-specific routines, if they are available !
	403	!----------------------------------------------------------------------!
	404
	405	IMPLICIT NONE
	406
	407	CHARACTER (LEN=*) :: direction
	408	INTEGER :: j, jshape(1)
	409
	410	!kk REAL, DIMENSION(:) :: ar !kk Does NOT work (Bug??)
	411	REAL, DIMENSION(0:ny) :: ar
	412	REAL, DIMENSION(0:ny+2) :: work
	413	REAL, DIMENSION(ny+2) :: work1
	414	COMPLEX, DIMENSION(:), ALLOCATABLE :: cwork
	415	#if defined( __ibm )
	416	REAL, DIMENSION(nau2) :: auy2, auy4
	417	#elif defined( __nec )
	418	REAL, DIMENSION(6*(ny+1)) :: work2
	419	#endif
	420
	421	IF ( fft_method == 'singleton-algorithm' ) THEN
	422
	423	!
	424	!-- Performing the fft with singleton's software works on every system,
	425	!-- since it is part of the model
	426	ALLOCATE( cwork(0:ny) )
	427
	428	IF ( direction == 'forward') THEN
	429
	430	DO j = 0, ny
	431	cwork(j) = CMPLX( ar(j) )
	432	ENDDO
	433
	434	jshape = SHAPE( cwork )
	435	CALL FFTN( cwork, jshape )
	436
	437	DO j = 0, (ny+1)/2
	438	ar(j) = REAL( cwork(j) )
	439	ENDDO
	440	DO j = 1, (ny+1)/2 - 1
	441	ar(ny+1-j) = -AIMAG( cwork(j) )
	442	ENDDO
	443
	444	ELSE
	445
	446	cwork(0) = CMPLX( ar(0), 0.0 )
	447	DO j = 1, (ny+1)/2 - 1
	448	cwork(j) = CMPLX( ar(j), -ar(ny+1-j) )
	449	cwork(ny+1-j) = CMPLX( ar(j), ar(ny+1-j) )
	450	ENDDO
	451	cwork((ny+1)/2) = CMPLX( ar((ny+1)/2), 0.0 )
	452
	453	jshape = SHAPE( cwork )
	454	CALL FFTN( cwork, jshape, inv = .TRUE. )
	455
	456	DO j = 0, ny
	457	ar(j) = REAL( cwork(j) )
	458	ENDDO
	459
	460	ENDIF
	461
	462	DEALLOCATE( cwork )
	463
	464	ELSEIF ( fft_method == 'temperton-algorithm' ) THEN
	465
	466	!
	467	!-- Performing the fft with Temperton's software works on every system,
	468	!-- since it is part of the model
	469	IF ( direction == 'forward' ) THEN
	470
	471	work(0:ny) = ar
	472	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
	473
	474	DO j = 0, (ny+1)/2
	475	ar(j) = work(2*j)
	476	ENDDO
	477	DO j = 1, (ny+1)/2 - 1
	478	ar(ny+1-j) = work(2*j+1)
	479	ENDDO
	480
	481	ELSE
	482
	483	DO j = 0, (ny+1)/2
	484	work(2*j) = ar(j)
	485	ENDDO
	486	DO j = 1, (ny+1)/2 - 1
	487	work(2*j+1) = ar(ny+1-j)
	488	ENDDO
	489	work(1) = 0.0
	490	work(ny+2) = 0.0
	491
	492	CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
	493	ar = work(0:ny)
	494
	495	ENDIF
	496
	497	ELSEIF ( fft_method == 'system-specific' ) THEN
	498
	499	#if defined( __ibm ) && ! defined( __ibmy_special )
	500	IF ( direction == 'forward') THEN
	501
	502	CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_ny, auy1, nau1, &
	503	auy2, nau2 )
	504
	505	DO j = 0, (ny+1)/2
	506	ar(j) = work(2*j)
	507	ENDDO
	508	DO j = 1, (ny+1)/2 - 1
	509	ar(ny+1-j) = work(2*j+1)
	510	ENDDO
	511
	512	ELSE
	513
	514	DO j = 0, (ny+1)/2
	515	work(2*j) = ar(j)
	516	ENDDO
	517	DO j = 1, (ny+1)/2 - 1
	518	work(2*j+1) = ar(ny+1-j)
	519	ENDDO
	520	work(1) = 0.0
	521	work(ny+2) = 0.0
	522
	523	CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_ny, auy3, nau1, &
	524	auy4, nau2 )
	525
	526	DO j = 0, ny
	527	ar(j) = work(j)
	528	ENDDO
	529
	530	ENDIF
	531	#elif defined( __nec )
	532	IF ( direction == 'forward' ) THEN
	533
	534	work(0:ny) = ar(0:ny)
	535
	536	CALL DZFFT( 1, ny+1, sqr_ny, work, work, trig_yf, work2, 0 )
	537
	538	DO j = 0, (ny+1)/2
	539	ar(j) = work(2*j)
	540	ENDDO
	541	DO j = 1, (ny+1)/2 - 1
	542	ar(ny+1-j) = work(2*j+1)
	543	ENDDO
	544
	545	ELSE
	546
	547	DO j = 0, (ny+1)/2
	548	work(2*j) = ar(j)
	549	ENDDO
	550	DO j = 1, (ny+1)/2 - 1
	551	work(2*j+1) = ar(ny+1-j)
	552	ENDDO
	553	work(1) = 0.0
	554	work(ny+2) = 0.0
	555
	556	CALL ZDFFT( -1, ny+1, sqr_ny, work, work, trig_yb, work2, 0 )
	557
	558	ar(0:ny) = work(0:ny)
	559
	560	ENDIF
	561	#else
[254]	562	message_string = 'no system-specific fft-call available'
	563	CALL message( 'fft_y', 'PA0188', 1, 2, 0, 6, 0 )
	564
[1]	565	#endif
	566
	567	ELSE
	568
[274]	569	message_string = 'fft method "' // TRIM( fft_method) // &
	570	'" not available'
[254]	571	CALL message( 'fft_y', 'PA0189', 1, 2, 0, 6, 0 )
[1]	572
	573	ENDIF
	574
	575	END SUBROUTINE fft_y
	576
	577	SUBROUTINE fft_x_m( ar, direction )
	578
	579	!----------------------------------------------------------------------!
	580	! fft_x_m !
	581	! !
	582	! Fourier-transformation along x-direction !
	583	! Version for 1d domain decomposition !
	584	! using multiple 1D FFT from Math Keisan on NEC !
	585	! or Temperton-algorithm !
	586	! (no singleton-algorithm on NEC because it does not vectorize) !
	587	! !
	588	!----------------------------------------------------------------------!
	589
	590	IMPLICIT NONE
	591
	592	CHARACTER (LEN=*) :: direction
	593	INTEGER :: i, k, siza, sizw
	594
	595	REAL, DIMENSION(0:nx,nz) :: ar
	596	REAL, DIMENSION(0:nx+3,nz+1) :: ai
	597	REAL, DIMENSION(6*(nx+4),nz+1) :: work1
	598	#if defined( __nec )
	599	COMPLEX, DIMENSION((nx+4)/2+1,nz+1) :: work
	600	#endif
	601
	602	IF ( fft_method == 'temperton-algorithm' ) THEN
	603
	604	siza = SIZE( ai, 1 )
	605
	606	IF ( direction == 'forward') THEN
	607
	608	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	609	ai(nx+1:,:) = 0.0
	610
	611	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, -1 )
	612
	613	DO k = 1, nz
	614	DO i = 0, (nx+1)/2
	615	ar(i,k) = ai(2*i,k)
	616	ENDDO
	617	DO i = 1, (nx+1)/2 - 1
	618	ar(nx+1-i,k) = ai(2*i+1,k)
	619	ENDDO
	620	ENDDO
	621
	622	ELSE
	623
	624	DO k = 1, nz
	625	DO i = 0, (nx+1)/2
	626	ai(2*i,k) = ar(i,k)
	627	ENDDO
	628	DO i = 1, (nx+1)/2 - 1
	629	ai(2*i+1,k) = ar(nx+1-i,k)
	630	ENDDO
	631	ai(1,k) = 0.0
	632	ai(nx+2,k) = 0.0
	633	ENDDO
	634
	635	CALL fft991cy( ai, work1, trigs_x, ifax_x, 1, siza, nx+1, nz, 1 )
	636
	637	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	638
	639	ENDIF
	640
	641	ELSEIF ( fft_method == 'system-specific' ) THEN
	642
	643	#if defined( __nec )
	644	siza = SIZE( ai, 1 )
	645	sizw = SIZE( work, 1 )
	646
	647	IF ( direction == 'forward') THEN
	648
	649	!
	650	!-- Tables are initialized once more. This call should not be
	651	!-- necessary, but otherwise program aborts in asymmetric case
	652	CALL DZFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
	653	trig_xf, work1, 0 )
	654
	655	ai(0:nx,1:nz) = ar(0:nx,1:nz)
	656	IF ( nz1 > nz ) THEN
	657	ai(:,nz1) = 0.0
	658	ENDIF
	659
	660	CALL DZFFTM( 1, nx+1, nz1, sqr_nx, ai, siza, work, sizw, &
	661	trig_xf, work1, 0 )
	662
	663	DO k = 1, nz
	664	DO i = 0, (nx+1)/2
	665	ar(i,k) = REAL( work(i+1,k) )
	666	ENDDO
	667	DO i = 1, (nx+1)/2 - 1
	668	ar(nx+1-i,k) = AIMAG( work(i+1,k) )
	669	ENDDO
	670	ENDDO
	671
	672	ELSE
	673
	674	!
	675	!-- Tables are initialized once more. This call should not be
	676	!-- necessary, but otherwise program aborts in asymmetric case
	677	CALL ZDFFTM( 0, nx+1, nz1, sqr_nx, work, nx+4, work, nx+4, &
	678	trig_xb, work1, 0 )
	679
	680	IF ( nz1 > nz ) THEN
	681	work(:,nz1) = 0.0
	682	ENDIF
	683	DO k = 1, nz
	684	work(1,k) = CMPLX( ar(0,k), 0.0 )
	685	DO i = 1, (nx+1)/2 - 1
	686	work(i+1,k) = CMPLX( ar(i,k), ar(nx+1-i,k) )
	687	ENDDO
	688	work(((nx+1)/2)+1,k) = CMPLX( ar((nx+1)/2,k), 0.0 )
	689	ENDDO
	690
	691	CALL ZDFFTM( -1, nx+1, nz1, sqr_nx, work, sizw, ai, siza, &
	692	trig_xb, work1, 0 )
	693
	694	ar(0:nx,1:nz) = ai(0:nx,1:nz)
	695
	696	ENDIF
	697
	698	#else
[254]	699	message_string = 'no system-specific fft-call available'
	700	CALL message( 'fft_x_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	701	#endif
	702
	703	ELSE
	704
[274]	705	message_string = 'fft method "' // TRIM( fft_method) // &
	706	'" not available'
[254]	707	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	708
	709	ENDIF
	710
	711	END SUBROUTINE fft_x_m
	712
	713	SUBROUTINE fft_y_m( ar, ny1, direction )
	714
	715	!----------------------------------------------------------------------!
	716	! fft_y_m !
	717	! !
	718	! Fourier-transformation along y-direction !
	719	! Version for 1d domain decomposition !
	720	! using multiple 1D FFT from Math Keisan on NEC !
	721	! or Temperton-algorithm !
	722	! (no singleton-algorithm on NEC because it does not vectorize) !
	723	! !
	724	!----------------------------------------------------------------------!
	725
	726	IMPLICIT NONE
	727
	728	CHARACTER (LEN=*) :: direction
	729	INTEGER :: j, k, ny1, siza, sizw
	730
	731	REAL, DIMENSION(0:ny1,nz) :: ar
	732	REAL, DIMENSION(0:ny+3,nz+1) :: ai
	733	REAL, DIMENSION(6*(ny+4),nz+1) :: work1
	734	#if defined( __nec )
	735	COMPLEX, DIMENSION((ny+4)/2+1,nz+1) :: work
	736	#endif
	737
	738	IF ( fft_method == 'temperton-algorithm' ) THEN
	739
	740	siza = SIZE( ai, 1 )
	741
	742	IF ( direction == 'forward') THEN
	743
	744	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	745	ai(ny+1:,:) = 0.0
	746
	747	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, -1 )
	748
	749	DO k = 1, nz
	750	DO j = 0, (ny+1)/2
	751	ar(j,k) = ai(2*j,k)
	752	ENDDO
	753	DO j = 1, (ny+1)/2 - 1
	754	ar(ny+1-j,k) = ai(2*j+1,k)
	755	ENDDO
	756	ENDDO
	757
	758	ELSE
	759
	760	DO k = 1, nz
	761	DO j = 0, (ny+1)/2
	762	ai(2*j,k) = ar(j,k)
	763	ENDDO
	764	DO j = 1, (ny+1)/2 - 1
	765	ai(2*j+1,k) = ar(ny+1-j,k)
	766	ENDDO
	767	ai(1,k) = 0.0
	768	ai(ny+2,k) = 0.0
	769	ENDDO
	770
	771	CALL fft991cy( ai, work1, trigs_y, ifax_y, 1, siza, ny+1, nz, 1 )
	772
	773	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	774
	775	ENDIF
	776
	777	ELSEIF ( fft_method == 'system-specific' ) THEN
	778
	779	#if defined( __nec )
	780	siza = SIZE( ai, 1 )
	781	sizw = SIZE( work, 1 )
	782
	783	IF ( direction == 'forward') THEN
	784
	785	!
	786	!-- Tables are initialized once more. This call should not be
	787	!-- necessary, but otherwise program aborts in asymmetric case
	788	CALL DZFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
	789	trig_yf, work1, 0 )
	790
	791	ai(0:ny,1:nz) = ar(0:ny,1:nz)
	792	IF ( nz1 > nz ) THEN
	793	ai(:,nz1) = 0.0
	794	ENDIF
	795
	796	CALL DZFFTM( 1, ny+1, nz1, sqr_ny, ai, siza, work, sizw, &
	797	trig_yf, work1, 0 )
	798
	799	DO k = 1, nz
	800	DO j = 0, (ny+1)/2
	801	ar(j,k) = REAL( work(j+1,k) )
	802	ENDDO
	803	DO j = 1, (ny+1)/2 - 1
	804	ar(ny+1-j,k) = AIMAG( work(j+1,k) )
	805	ENDDO
	806	ENDDO
	807
	808	ELSE
	809
	810	!
	811	!-- Tables are initialized once more. This call should not be
	812	!-- necessary, but otherwise program aborts in asymmetric case
	813	CALL ZDFFTM( 0, ny+1, nz1, sqr_ny, work, ny+4, work, ny+4, &
	814	trig_yb, work1, 0 )
	815
	816	IF ( nz1 > nz ) THEN
	817	work(:,nz1) = 0.0
	818	ENDIF
	819	DO k = 1, nz
	820	work(1,k) = CMPLX( ar(0,k), 0.0 )
	821	DO j = 1, (ny+1)/2 - 1
	822	work(j+1,k) = CMPLX( ar(j,k), ar(ny+1-j,k) )
	823	ENDDO
	824	work(((ny+1)/2)+1,k) = CMPLX( ar((ny+1)/2,k), 0.0 )
	825	ENDDO
	826
	827	CALL ZDFFTM( -1, ny+1, nz1, sqr_ny, work, sizw, ai, siza, &
	828	trig_yb, work1, 0 )
	829
	830	ar(0:ny,1:nz) = ai(0:ny,1:nz)
	831
	832	ENDIF
	833
	834	#else
[254]	835	message_string = 'no system-specific fft-call available'
	836	CALL message( 'fft_y_m', 'PA0188', 1, 2, 0, 6, 0 )
[1]	837	#endif
	838
	839	ELSE
[254]	840
[274]	841	message_string = 'fft method "' // TRIM( fft_method) // &
	842	'" not available'
[254]	843	CALL message( 'fft_x_m', 'PA0189', 1, 2, 0, 6, 0 )
[1]	844
	845	ENDIF
	846
	847	END SUBROUTINE fft_y_m
	848
	849	END MODULE fft_xy

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |