Home

Context Navigation

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

Last change on this file since 357 was 274, checked in by heinze, 16 years ago
Indentation of the message calls corrected
Property svn:keywords set to `Id`
File size: 24.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |