Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |