Home

Context Navigation

fft_xy.f90 @ 961

Last change on this file since 961 was 392, checked in by raasch, 15 years ago

New:
---

Adapted for machine lck
(mrun, mbuild, subjob)

bc_lr/bc_ns in most subroutines replaced by LOGICAL variables bc_lr_cyc,
bc_ns_cyc for speed optimization
(check_parameters, diffusion_u, diffusion_v, diffusion_w, modules)

Additional timestep criterion in case of simulations with plant canopy (timestep)

Check for illegal entries in section_xy|xz|yz that exceed nz+1|ny+1|nx+1
(check_parameters)

Clipping of dvrp output implemented. Default colourtable for particles
implemented, particle attributes (color, dvrp_size) can be set with new
parameters particle_color, particle_dvrpsize, color_interval,
dvrpsize_interval (init_dvrp, data_output_dvrp, modules, user_data_output_dvrp).
Slicer attributes (dvrp) are set with new routine set_slicer_attributes_dvrp
and are controlled with existing parameters slicer_range_limits.
(set_slicer_attributes_dvrp)

Ocean atmosphere coupling allows to use independent precursor runs in order
to account for different spin-up times. The time when coupling has to be
started is given by new inipar parameter coupling_start_time. The precursor
ocean run has to be started using new mrun option "-y" in order to add
appendix "_O" to all output files.
(check_for_restart, check_parameters, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, header, init_coupling, modules, mrun,
parin, read_var_list, surface_coupler, time_integration, write_var_list)

Polygon reduction for topography and ground plate isosurface. Reduction level
for buildings can be chosen with parameter cluster_size. (init_dvrp)

External pressure gradient (check_parameters, header, init_3d_model, modules,
parin, prognostic_equations, read_var_list, write_var_list)

New topography case 'single_street_canyon' (header, init_grid, modules, parin,
read_var_list, user_check_parameters, user_header, user_init_grid, write_var_list)

Option to predefine a target bulk velocity for conserve_volume_flow
(check_parameters, header, init_3d_model, modules, parin, read_var_list,
write_var_list)

Option for user defined 2D data output in xy cross sections at z=nzb+1
(data_output_2d, user_data_output_2d)

xy cross section output of surface heatfluxes (latent, sensible)
(average_3d_data, check_parameters, data_output_2d, modules, read_3d_binary,
sum_up_3d_data, write_3d_binary)

average_3d_data, check_for_restart, check_parameters, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, init_coupling, init_dvrp, init_grid, init_3d_model, header, mbuild, modules, mrun, package_parin, parin, prognostic_equations, read_3d_binary, read_var_list, subjob, surface_coupler, timestep, time_integration, user_check_parameters, user_data_output_2d, user_data_output_dvrp, user_header, user_init_grid, write_3d_binary, write_var_list

New: set_particle_attributes, set_slicer_attributes_dvrp

Changed:

lcmuk changed to lc to avoid problems with Intel compiler on sgi-ice
(poisfft)

For extended NetCDF files, the updated title attribute includes an update of
time_average_text where appropriate. (netcdf)

In case of restart runs without extension, initial profiles are not written
to NetCDF-file anymore. (data_output_profiles, modules, read_var_list, write_var_list)

Small change in formatting of the message handling routine concering the output in the
job protocoll. (message)

initializing_actions='read_data_for_recycling' renamed to 'cyclic_fill', now
independent of turbulent_inflow (check_parameters, header, init_3d_model)

2 NetCDF error numbers changed. (data_output_3d)

A Link to the website appendix_a.html is printed for further information
about the possible errors. (message)

Temperature gradient criterion for estimating the boundary layer height
replaced by the gradient criterion of Sullivan et al. (1998). (flow_statistics)

NetCDF unit attribute in timeseries output in case of statistic regions added
(netcdf)

Output of NetCDF messages with aid of message handling routine.
(check_open, close_file, data_output_2d, data_output_3d,
data_output_profiles, data_output_ptseries, data_output_spectra,
data_output_tseries, netcdf, output_particles_netcdf)

Output of messages replaced by message handling routine.
(advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart,
check_open, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp,
data_output_profiles, data_output_spectra, fft_xy, flow_statistics, header,
init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid,
netcdf, parin, plant_canopy_model, poisfft_hybrid, poismg, read_3d_binary,
read_var_list, surface_coupler, temperton_fft, timestep, user_actions,
user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy,
user_parin, user_read_restart_data, user_spectra )

Maximum number of tails is calculated from maximum number of particles and
skip_particles_for_tail (init_particles)

Value of vertical_particle_advection may differ for each particle group
(advec_particles, header, modules)

First constant in array den also defined as type double. (eqn_state_seawater)

Parameter dvrp_psize moved from particles_par to dvrp_graphics_par. (package_parin)

topography_grid_convention moved from userpar to inipar (check_parameters,
header, parin, read_var_list, user_check_parameters, user_header,
user_init_grid, user_parin, write_var_list)

Default value of grid_matching changed to strict.

Adjustments for runs on lcxt4 (necessary due to an software update on CRAY) and
for coupled runs on ibmy (mrun, subjob)

advec_particles, advec_s_bc, buoyancy, calc_spectra, check_for_restart, check_open, check_parameters, close_file, coriolis, cpu_log, data_output_2d, data_output_3d, data_output_dvrp, data_output_profiles, data_output_ptseries, data_output_spectra, data_output_tseries, eqn_state_seawater, fft_xy, flow_statistics, header, init_1d_model, init_3d_model, init_dvrp, init_grid, init_particles, init_pegrid, message, mrun, netcdf, output_particles_netcdf, package_parin, parin, plant_canopy_model, poisfft, poisfft_hybrid, poismg, read_3d_binary, read_var_list, sort_particles, subjob, user_check_parameters, user_header, user_init_grid, user_parin, surface_coupler, temperton_fft, timestep, user_actions, user_data_output_dvrp, user_dvrp_coltab, user_init_grid, user_init_plant_canopy, user_parin, user_read_restart_data, user_spectra, write_var_list

Errors:

Bugfix: Initial hydrostatic pressure profile in case of ocean runs is now
calculated in 5 iteration steps. (init_ocean)

Bugfix: wrong sign in buoyancy production of ocean part in case of not using
the reference density (only in 3D routine production_e) (production_e)

Bugfix: output of averaged 2d/3d quantities requires that an avaraging
interval has been set, respective error message is included (check_parameters)

Bugfix: Output on unit 14 only if requested by write_binary.
(user_last_actions)

Bugfix to avoid zero division by km_neutral (production_e)

Bugfix for extended NetCDF files: In order to avoid 'data mode' errors if
updated attributes are larger than their original size, NF90_PUT_ATT is called
in 'define mode' enclosed by NF90_REDEF and NF90_ENDDEF calls. This implies a
possible performance loss; an alternative strategy would be to ensure equal
attribute size in a job chain. (netcdf)

Bugfix: correction of initial volume flow for non-flat topography (init_3d_model)
Bugfix: zero initialization of arrays within buildings for 'cyclic_fill' (init_3d_model)

Bugfix: to_be_resorted => s_av for time-averaged scalars (data_output_2d, data_output_3d)

Bugfix: error in formatting the output (message)

Bugfix: avoid that ngp_2dh_s_inner becomes zero (init_3_model)

Typographical error: unit of wpt in dots_unit (modules)

Bugfix: error in check, if particles moved further than one subdomain length.
This check must not be applied for newly released particles. (advec_particles)

Bugfix: several tail counters are initialized, particle_tail_coordinates is
only written to file if its third index is > 0, arrays for tails are allocated
with a minimum size of 10 tails if there is no tail initially (init_particles,
advec_particles)

Bugfix: pressure included for profile output (check_parameters)

Bugfix: Type of count and count_rate changed to default INTEGER on NEC machines
(cpu_log)

Bugfix: output if particle time series only if particle advection is switched
on. (time_integration)

Bugfix: qsws was calculated in case of constant heatflux = .FALSE. (prandtl_fluxes)

Bugfix: averaging along z is not allowed for 2d quantities (e.g. u* and z0) (data_output_2d)

Typographical errors (netcdf)

If the inversion height calculated by the prerun is zero, inflow_damping_height
must be explicitly specified (init_3d_model)

Small bugfix concerning 3d 64bit netcdf output format (header)

Bugfix: dt_fixed removed from the restart file, because otherwise, no change
from a fixed to a variable timestep would be possible in restart runs.
(read_var_list, write_var_list)

Bugfix: initial setting of time_coupling in coupled restart runs (time_integration)

advec_particles, check_parameters, cpu_log, data_output_2d, data_output_3d, header, init_3d_model, init_particles, init_ocean, modules, netcdf, prandtl_fluxes, production_e, read_var_list, time_integration, user_last_actions, write_var_list

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

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

Context Navigation

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

Download in other formats: