Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |