source: palm/trunk/SOURCE/poisfft.f90 @ 392

Last change on this file since 392 was 392, checked in by raasch, 12 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: 44.6 KB
Line 
1 MODULE poisfft_mod
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: poisfft.f90 392 2009-09-24 10:39:14Z raasch $
11!
12! 377 2009-09-04 11:09:00Z raasch
13! __lcmuk changed to __lc to avoid problems with Intel compiler on sgi-ice
14!
15! 164 2008-05-15 08:46:15Z raasch
16! Arguments removed from transpose routines
17!
18! 128 2007-10-26 13:11:14Z raasch
19! Bugfix: wavenumber calculation for even nx in routines maketri
20!
21! 85 2007-05-11 09:35:14Z raasch
22! Bugfix: work_fft*_vec removed from some PRIVATE-declarations
23!
24! 76 2007-03-29 00:58:32Z raasch
25! Tridiagonal coefficients adjusted for Neumann boundary conditions both at
26! the bottom and the top.
27!
28! RCS Log replace by Id keyword, revision history cleaned up
29!
30! Revision 1.24  2006/08/04 15:00:24  raasch
31! Default setting of the thread number tn in case of not using OpenMP
32!
33! Revision 1.23  2006/02/23 12:48:38  raasch
34! Additional compiler directive in routine tridia_1dd for preventing loop
35! exchange on NEC-SX6
36!
37! Revision 1.20  2004/04/30 12:38:09  raasch
38! Parts of former poisfft_hybrid moved to this subroutine,
39! former subroutine changed to a module, renaming of FFT-subroutines and
40! -module, FFTs completely substituted by calls of fft_x and fft_y,
41! NAG fft used in the non-parallel case completely removed, l in maketri
42! is now a 1d-array, variables passed by modules instead of using parameter
43! lists, enlarged transposition arrays introduced
44!
45! Revision 1.1  1997/07/24 11:24:14  raasch
46! Initial revision
47!
48!
49! Description:
50! ------------
51! See below.
52!------------------------------------------------------------------------------!
53
54!--------------------------------------------------------------------------!
55!                             poisfft                                      !
56!                                                                          !
57!                Original version: Stephan Siano (pois3d)                  !
58!                                                                          !
59!  Institute of Meteorology and Climatology, University of Hannover        !
60!                             Germany                                      !
61!                                                                          !
62!  Version as of July 23,1996                                              !
63!                                                                          !
64!                                                                          !
65!        Version for parallel computers: Siegfried Raasch                  !
66!                                                                          !
67!  Version as of July 03,1997                                              !
68!                                                                          !
69!  Solves the Poisson equation with a 2D spectral method                   !
70!        d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s                    !
71!                                                                          !
72!  Input:                                                                  !
73!  real    ar                 contains in the (nnx,nny,nnz) elements,      !
74!                             starting from the element (1,nys,nxl), the   !
75!                             values for s                                 !
76!  real    work               Temporary array                              !
77!                                                                          !
78!  Output:                                                                 !
79!  real    ar                 contains the solution for p                  !
80!--------------------------------------------------------------------------!
81
82    USE fft_xy
83    USE indices
84    USE transpose_indices
85
86    IMPLICIT NONE
87
88    PRIVATE
89    PUBLIC  poisfft, poisfft_init
90
91    INTERFACE poisfft
92       MODULE PROCEDURE poisfft
93    END INTERFACE poisfft
94
95    INTERFACE poisfft_init
96       MODULE PROCEDURE poisfft_init
97    END INTERFACE poisfft_init
98
99 CONTAINS
100
101    SUBROUTINE poisfft_init
102
103       CALL fft_init
104
105    END SUBROUTINE poisfft_init
106
107
108    SUBROUTINE poisfft( ar, work )
109
110       USE cpulog
111       USE interfaces
112       USE pegrid
113
114       IMPLICIT NONE
115
116       REAL, DIMENSION(1:nza,nys:nyna,nxl:nxra) ::  ar, work
117
118
119       CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
120
121!
122!--    Two-dimensional Fourier Transformation in x- and y-direction.
123#if defined( __parallel )
124       IF ( pdims(2) == 1 )  THEN
125
126!
127!--       1d-domain-decomposition along x:
128!--       FFT along y and transposition y --> x
129          CALL ffty_tr_yx( ar, work, ar )
130
131!
132!--       FFT along x, solving the tridiagonal system and backward FFT
133          CALL fftx_tri_fftx( ar )
134
135!
136!--       Transposition x --> y and backward FFT along y
137          CALL tr_xy_ffty( ar, work, ar )
138
139       ELSEIF ( pdims(1) == 1 )  THEN
140
141!
142!--       1d-domain-decomposition along y:
143!--       FFT along x and transposition x --> y
144          CALL fftx_tr_xy( ar, work, ar )
145
146!
147!--       FFT along y, solving the tridiagonal system and backward FFT
148          CALL ffty_tri_ffty( ar )
149
150!
151!--       Transposition y --> x and backward FFT along x
152          CALL tr_yx_fftx( ar, work, ar )
153
154       ELSE
155
156!
157!--       2d-domain-decomposition
158!--       Transposition z --> x
159          CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
160          CALL transpose_zx( ar, work, ar )
161          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
162
163          CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
164          CALL fftxp( ar, 'forward' )
165          CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
166
167!
168!--       Transposition x --> y
169          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
170          CALL transpose_xy( ar, work, ar )
171          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
172
173          CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
174          CALL fftyp( ar, 'forward' )
175          CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
176
177!
178!--       Transposition y --> z
179          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
180          CALL transpose_yz( ar, work, ar )
181          CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
182
183!
184!--       Solve the Poisson equation in z-direction in cartesian space.
185          CALL cpu_log( log_point_s(6), 'tridia', 'start' )
186          CALL tridia( ar )
187          CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
188
189!
190!--       Inverse Fourier Transformation
191!--       Transposition z --> y
192          CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
193          CALL transpose_zy( ar, work, ar )
194          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
195
196          CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
197          CALL fftyp( ar, 'backward' )
198          CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
199
200!
201!--       Transposition y --> x
202          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
203          CALL transpose_yx( ar, work, ar )
204          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
205
206          CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
207          CALL fftxp( ar, 'backward' )
208          CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
209
210!
211!--       Transposition x --> z
212          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
213          CALL transpose_xz( ar, work, ar )
214          CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
215
216       ENDIF
217
218#else
219
220!
221!--    Two-dimensional Fourier Transformation along x- and y-direction.
222       CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
223       CALL fftx( ar, 'forward' )
224       CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
225       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
226       CALL ffty( ar, 'forward' )
227       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
228
229!
230!--    Solve the Poisson equation in z-direction in cartesian space.
231       CALL cpu_log( log_point_s(6), 'tridia', 'start' )
232       CALL tridia( ar )
233       CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
234
235!
236!--    Inverse Fourier Transformation.
237       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
238       CALL ffty( ar, 'backward' )
239       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
240       CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
241       CALL fftx( ar, 'backward' )
242       CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
243
244#endif
245
246       CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
247
248    END SUBROUTINE poisfft
249
250
251
252    SUBROUTINE tridia( ar )
253
254!------------------------------------------------------------------------------!
255! solves the linear system of equations:
256!
257! -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
258!   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
259! 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
260!
261! by using the Thomas algorithm
262!------------------------------------------------------------------------------!
263
264       USE arrays_3d
265
266       IMPLICIT NONE
267
268       INTEGER ::  i, j, k, nnyh
269
270       REAL, DIMENSION(nxl_z:nxr_z,0:nz-1)   ::  ar1
271       REAL, DIMENSION(5,nxl_z:nxr_z,0:nz-1) ::  tri
272
273#if defined( __parallel )
274       REAL    ::  ar(nxl_z:nxr_za,nys_z:nyn_za,1:nza)
275#else
276       REAL    ::  ar(1:nz,nys_z:nyn_z,nxl_z:nxr_z)
277#endif
278
279
280       nnyh = (ny+1) / 2
281
282!
283!--    Define constant elements of the tridiagonal matrix.
284       DO  k = 0, nz-1
285          DO  i = nxl_z, nxr_z
286             tri(2,i,k) = ddzu(k+1) * ddzw(k+1)
287             tri(3,i,k) = ddzu(k+2) * ddzw(k+1)
288          ENDDO
289       ENDDO
290
291#if defined( __parallel )
292!
293!--    Repeat for all y-levels.
294       DO  j = nys_z, nyn_z
295          IF ( j <= nnyh )  THEN
296             CALL maketri( tri, j )
297          ELSE
298             CALL maketri( tri, ny+1-j )
299          ENDIF
300          CALL split( tri )
301          CALL substi( ar, ar1, tri, j )
302       ENDDO
303#else
304!
305!--    First y-level.
306       CALL maketri( tri, nys_z )
307       CALL split( tri )
308       CALL substi( ar, ar1, tri, 0 )
309
310!
311!--    Further y-levels.
312       DO  j = 1, nnyh - 1
313          CALL maketri( tri, j )
314          CALL split( tri )
315          CALL substi( ar, ar1, tri, j )
316          CALL substi( ar, ar1, tri, ny+1-j )
317       ENDDO
318       CALL maketri( tri, nnyh )
319       CALL split( tri )
320       CALL substi( ar, ar1, tri, nnyh+nys )
321#endif
322
323    CONTAINS
324
325       SUBROUTINE maketri( tri, j )
326
327!------------------------------------------------------------------------------!
328! Computes the i- and j-dependent component of the matrix
329!------------------------------------------------------------------------------!
330
331          USE arrays_3d
332          USE constants
333          USE control_parameters
334          USE grid_variables
335
336          IMPLICIT NONE
337
338          INTEGER ::  i, j, k, nnxh
339          REAL    ::  a, c
340          REAL    ::  ll(nxl_z:nxr_z)
341          REAL    ::  tri(5,nxl_z:nxr_z,0:nz-1)
342
343
344          nnxh = ( nx + 1 ) / 2
345
346!
347!--       Provide the tridiagonal matrix for solution of the Poisson equation in
348!--       Fourier space. The coefficients are computed following the method of
349!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
350!--       Siano's original version by discretizing the Poisson equation,
351!--       before it is Fourier-transformed
352#if defined( __parallel )
353          DO  i = nxl_z, nxr_z
354             IF ( i >= 0 .AND. i <= nnxh )  THEN
355                ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
356                                          FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
357                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
358                                          FLOAT( ny+1 ) ) ) / ( dy * dy )
359             ELSE
360                ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
361                                          FLOAT( nx+1 ) ) ) / ( dx * dx ) + &
362                        2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
363                                          FLOAT( ny+1 ) ) ) / ( dy * dy )
364             ENDIF
365             DO  k = 0,nz-1
366                a = -1.0 * ddzu(k+2) * ddzw(k+1)
367                c = -1.0 * ddzu(k+1) * ddzw(k+1)
368                tri(1,i,k) = a + c - ll(i)
369             ENDDO
370          ENDDO
371#else
372          DO  i = 0, nnxh
373             ll(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / FLOAT( nx+1 ) ) ) / &
374                           ( dx * dx ) + &
375                     2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / FLOAT( ny+1 ) ) ) / &
376                           ( dy * dy )
377             DO  k = 0, nz-1
378                a = -1.0 * ddzu(k+2) * ddzw(k+1)
379                c = -1.0 * ddzu(k+1) * ddzw(k+1)
380                tri(1,i,k) = a + c - ll(i)
381                IF ( i >= 1 .and. i < nnxh )  THEN
382                   tri(1,nx+1-i,k) = tri(1,i,k)
383                ENDIF
384             ENDDO
385          ENDDO
386#endif
387          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
388             DO  i = nxl_z, nxr_z
389                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
390             ENDDO
391          ENDIF
392          IF ( ibc_p_t == 1 )  THEN
393             DO  i = nxl_z, nxr_z
394                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
395             ENDDO
396          ENDIF
397
398       END SUBROUTINE maketri
399
400
401       SUBROUTINE substi( ar, ar1, tri, j )
402
403!------------------------------------------------------------------------------!
404! Substitution (Forward and Backward) (Thomas algorithm)
405!------------------------------------------------------------------------------!
406
407          USE control_parameters
408
409          IMPLICIT NONE
410
411          INTEGER ::  i, j, k
412          REAL    ::  ar1(nxl_z:nxr_z,0:nz-1)
413          REAL    ::  tri(5,nxl_z:nxr_z,0:nz-1)
414#if defined( __parallel )
415          REAL    ::  ar(nxl_z:nxr_za,nys_z:nyn_za,1:nza)
416#else
417          REAL    ::  ar(1:nz,nys_z:nyn_z,nxl_z:nxr_z)
418#endif
419
420!
421!--       Forward substitution.
422          DO  i = nxl_z, nxr_z
423#if defined( __parallel )
424             ar1(i,0) = ar(i,j,1)
425#else
426             ar1(i,0) = ar(1,j,i)
427#endif
428          ENDDO
429          DO  k = 1, nz - 1
430             DO  i = nxl_z, nxr_z
431#if defined( __parallel )
432                ar1(i,k) = ar(i,j,k+1) - tri(5,i,k) * ar1(i,k-1)
433#else
434                ar1(i,k) = ar(k+1,j,i) - tri(5,i,k) * ar1(i,k-1)
435#endif
436             ENDDO
437          ENDDO
438
439!
440!--       Backward substitution.
441          DO  i = nxl_z, nxr_z
442#if defined( __parallel )
443             ar(i,j,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
444#else
445             ar(nz,j,i) = ar1(i,nz-1) / tri(4,i,nz-1)
446#endif
447          ENDDO
448          DO  k = nz-2, 0, -1
449             DO  i = nxl_z, nxr_z
450#if defined( __parallel )
451                ar(i,j,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,j,k+2) ) &
452                              / tri(4,i,k)
453#else
454                ar(k+1,j,i) = ( ar1(i,k) - tri(3,i,k) * ar(k+2,j,i) ) &
455                              / tri(4,i,k)
456#endif
457             ENDDO
458          ENDDO
459
460!
461!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
462!--       The respective values of ar should be zero at all k-levels if
463!--       acceleration of horizontally averaged vertical velocity is zero.
464          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
465             IF ( j == 0  .AND.  nxl_z == 0 )  THEN
466#if defined( __parallel )
467                DO  k = 1, nz
468                   ar(nxl_z,j,k) = 0.0
469                ENDDO
470#else
471                DO  k = 1, nz
472                   ar(k,j,nxl_z) = 0.0
473                ENDDO
474#endif
475             ENDIF
476          ENDIF
477
478       END SUBROUTINE substi
479
480
481       SUBROUTINE split( tri )
482
483!------------------------------------------------------------------------------!
484! Splitting of the tridiagonal matrix (Thomas algorithm)
485!------------------------------------------------------------------------------!
486
487          IMPLICIT NONE
488
489          INTEGER ::  i, k
490          REAL    ::  tri(5,nxl_z:nxr_z,0:nz-1)
491
492!
493!--       Splitting.
494          DO  i = nxl_z, nxr_z
495             tri(4,i,0) = tri(1,i,0)
496          ENDDO
497          DO  k = 1, nz-1
498             DO  i = nxl_z, nxr_z
499                tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
500                tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
501             ENDDO
502          ENDDO
503
504       END SUBROUTINE split
505
506    END SUBROUTINE tridia
507
508
509#if defined( __parallel )
510    SUBROUTINE fftxp( ar, direction )
511
512!------------------------------------------------------------------------------!
513! Fourier-transformation along x-direction                 Parallelized version
514!------------------------------------------------------------------------------!
515
516       IMPLICIT NONE
517
518       CHARACTER (LEN=*) ::  direction
519       INTEGER           ::  j, k
520       REAL              ::  ar(0:nxa,nys_x:nyn_xa,nzb_x:nzt_xa)
521
522!
523!--    Performing the fft with one of the methods implemented
524       DO  k = nzb_x, nzt_x
525          DO  j = nys_x, nyn_x
526             CALL fft_x( ar(0:nx,j,k), direction )
527          ENDDO
528       ENDDO
529
530    END SUBROUTINE fftxp
531
532#else
533    SUBROUTINE fftx( ar, direction )
534
535!------------------------------------------------------------------------------!
536! Fourier-transformation along x-direction                 Non parallel version
537!------------------------------------------------------------------------------!
538
539       IMPLICIT NONE
540
541       CHARACTER (LEN=*) ::  direction
542       INTEGER           ::  i, j, k
543       REAL              ::  ar(1:nz,0:ny,0:nx)
544
545!
546!--    Performing the fft with one of the methods implemented
547       DO  k = 1, nz
548          DO  j = 0, ny
549             CALL fft_x( ar(k,j,0:nx), direction )
550          ENDDO
551       ENDDO
552
553    END SUBROUTINE fftx
554#endif
555
556
557#if defined( __parallel )
558    SUBROUTINE fftyp( ar, direction )
559
560!------------------------------------------------------------------------------!
561! Fourier-transformation along y-direction                 Parallelized version
562!------------------------------------------------------------------------------!
563
564       IMPLICIT NONE
565
566       CHARACTER (LEN=*) ::  direction
567       INTEGER           ::  i, k
568       REAL              ::  ar(0:nya,nxl_y:nxr_ya,nzb_y:nzt_ya)
569
570!
571!--    Performing the fft with one of the methods implemented
572       DO  k = nzb_y, nzt_y
573          DO  i = nxl_y, nxr_y
574             CALL fft_y( ar(0:ny,i,k), direction )
575          ENDDO
576       ENDDO
577
578    END SUBROUTINE fftyp
579
580#else
581    SUBROUTINE ffty( ar, direction )
582
583!------------------------------------------------------------------------------!
584! Fourier-transformation along y-direction                 Non parallel version
585!------------------------------------------------------------------------------!
586
587       IMPLICIT NONE
588
589       CHARACTER (LEN=*) ::  direction
590       INTEGER           ::  i, k
591       REAL              ::  ar(1:nz,0:ny,0:nx)
592
593!
594!--    Performing the fft with one of the methods implemented
595       DO  k = 1, nz
596          DO  i = 0, nx
597             CALL fft_y( ar(k,0:ny,i), direction )
598          ENDDO
599       ENDDO
600
601    END SUBROUTINE ffty
602#endif
603
604#if defined( __parallel )
605    SUBROUTINE ffty_tr_yx( f_in, work, f_out )
606
607!------------------------------------------------------------------------------!
608!  Fourier-transformation along y with subsequent transposition y --> x for
609!  a 1d-decomposition along x
610!
611!  ATTENTION: The performance of this routine is much faster on the NEC-SX6,
612!             if the first index of work_ffty_vec is odd. Otherwise
613!             memory bank conflicts may occur (especially if the index is a
614!             multiple of 128). That's why work_ffty_vec is dimensioned as
615!             0:ny+1.
616!             Of course, this will not work if users are using an odd number
617!             of gridpoints along y.
618!------------------------------------------------------------------------------!
619
620       USE control_parameters
621       USE cpulog
622       USE indices
623       USE interfaces
624       USE pegrid
625       USE transpose_indices
626
627       IMPLICIT NONE
628
629       INTEGER            ::  i, iend, iouter, ir, j, k
630       INTEGER, PARAMETER ::  stridex = 4
631
632       REAL, DIMENSION(0:ny,stridex)                    ::  work_ffty
633#if defined( __nec )
634       REAL, DIMENSION(0:ny+1,1:nz,nxl:nxr)             ::  work_ffty_vec
635#endif
636       REAL, DIMENSION(1:nza,0:nya,nxl:nxra)            ::  f_in
637       REAL, DIMENSION(nnx,1:nza,nys_x:nyn_xa,pdims(1)) ::  f_out
638       REAL, DIMENSION(nxl:nxra,1:nza,0:nya)            ::  work
639
640!
641!--    Carry out the FFT along y, where all data are present due to the
642!--    1d-decomposition along x. Resort the data in a way that x becomes
643!--    the first index.
644       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
645
646       IF ( host(1:3) == 'nec' )  THEN
647#if defined( __nec )
648!
649!--       Code optimized for vector processors
650!$OMP     PARALLEL PRIVATE ( i, j, k )
651!$OMP     DO
652          DO  i = nxl, nxr
653
654             DO  j = 0, ny
655                DO  k = 1, nz
656                   work_ffty_vec(j,k,i) = f_in(k,j,i)
657                ENDDO
658             ENDDO
659
660             CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
661
662          ENDDO
663
664!$OMP     DO
665          DO  k = 1, nz
666             DO  j = 0, ny
667                DO  i = nxl, nxr
668                   work(i,k,j) = work_ffty_vec(j,k,i)
669                ENDDO
670             ENDDO
671          ENDDO
672!$OMP     END PARALLEL
673#endif
674
675       ELSE
676
677!
678!--       Cache optimized code.
679!--       The i-(x-)direction is split into a strided outer loop and an inner
680!--       loop for better cache performance
681!$OMP     PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
682!$OMP     DO
683          DO  iouter = nxl, nxr, stridex
684
685             iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
686
687             DO  k = 1, nz
688
689                DO  i = iouter, iend
690
691                   ir = i-iouter+1  ! counter within a stride
692                   DO  j = 0, ny
693                      work_ffty(j,ir) = f_in(k,j,i)
694                   ENDDO
695!
696!--                FFT along y
697                   CALL fft_y( work_ffty(:,ir), 'forward' )
698
699                ENDDO
700
701!
702!--             Resort
703                DO  j = 0, ny
704                   DO  i = iouter, iend
705                      work(i,k,j) = work_ffty(j,i-iouter+1)
706                   ENDDO
707                ENDDO
708
709             ENDDO
710
711          ENDDO
712!$OMP     END PARALLEL
713
714       ENDIF
715       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
716
717!
718!--    Transpose array
719       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
720       CALL MPI_ALLTOALL( work(nxl,1,0),      sendrecvcount_xy, MPI_REAL, &
721                          f_out(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
722                          comm1dx, ierr )
723       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
724
725    END SUBROUTINE ffty_tr_yx
726
727
728    SUBROUTINE tr_xy_ffty( f_in, work, f_out )
729
730!------------------------------------------------------------------------------!
731!  Transposition x --> y with a subsequent backward Fourier transformation for
732!  a 1d-decomposition along x
733!------------------------------------------------------------------------------!
734
735       USE control_parameters
736       USE cpulog
737       USE indices
738       USE interfaces
739       USE pegrid
740       USE transpose_indices
741
742       IMPLICIT NONE
743
744       INTEGER            ::  i, iend, iouter, ir, j, k
745       INTEGER, PARAMETER ::  stridex = 4
746
747       REAL, DIMENSION(0:ny,stridex)                    ::  work_ffty
748#if defined( __nec )
749       REAL, DIMENSION(0:ny+1,1:nz,nxl:nxr)             ::  work_ffty_vec
750#endif
751       REAL, DIMENSION(nnx,1:nza,nys_x:nyn_xa,pdims(1)) ::  f_in
752       REAL, DIMENSION(1:nza,0:nya,nxl:nxra)            ::  f_out
753       REAL, DIMENSION(nxl:nxra,1:nza,0:nya)            ::  work
754
755!
756!--    Transpose array
757       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
758       CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, &
759                          work(nxl,1,0),     sendrecvcount_xy, MPI_REAL, &
760                          comm1dx, ierr )
761       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
762
763!
764!--    Resort the data in a way that y becomes the first index and carry out the
765!--    backward fft along y.
766       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
767
768       IF ( host(1:3) == 'nec' )  THEN
769#if defined( __nec )
770!
771!--       Code optimized for vector processors
772!$OMP     PARALLEL PRIVATE ( i, j, k )
773!$OMP     DO
774          DO  k = 1, nz
775             DO  j = 0, ny
776                DO  i = nxl, nxr
777                   work_ffty_vec(j,k,i) = work(i,k,j)
778                ENDDO
779             ENDDO
780          ENDDO
781
782!$OMP     DO
783          DO  i = nxl, nxr
784
785             CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
786
787             DO  j = 0, ny
788                DO  k = 1, nz
789                   f_out(k,j,i) = work_ffty_vec(j,k,i)
790                ENDDO
791             ENDDO
792
793          ENDDO
794!$OMP     END PARALLEL
795#endif
796
797       ELSE
798
799!
800!--       Cache optimized code.
801!--       The i-(x-)direction is split into a strided outer loop and an inner
802!--       loop for better cache performance
803!$OMP     PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
804!$OMP     DO
805          DO  iouter = nxl, nxr, stridex
806
807             iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
808
809             DO  k = 1, nz
810!
811!--             Resort
812                DO  j = 0, ny
813                   DO  i = iouter, iend
814                      work_ffty(j,i-iouter+1) = work(i,k,j)
815                   ENDDO
816                ENDDO
817
818                DO  i = iouter, iend
819
820!
821!--                FFT along y
822                   ir = i-iouter+1  ! counter within a stride
823                   CALL fft_y( work_ffty(:,ir), 'backward' )
824
825                   DO  j = 0, ny
826                      f_out(k,j,i) = work_ffty(j,ir)
827                   ENDDO
828                ENDDO
829
830             ENDDO
831
832          ENDDO
833!$OMP     END PARALLEL
834
835       ENDIF
836
837       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
838
839    END SUBROUTINE tr_xy_ffty
840
841
842    SUBROUTINE fftx_tri_fftx( ar )
843
844!------------------------------------------------------------------------------!
845!  FFT along x, solution of the tridiagonal system and backward FFT for
846!  a 1d-decomposition along x
847!
848!  WARNING: this subroutine may still not work for hybrid parallelization
849!           with OpenMP (for possible necessary changes see the original
850!           routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
851!------------------------------------------------------------------------------!
852
853       USE control_parameters
854       USE cpulog
855       USE grid_variables
856       USE indices
857       USE interfaces
858       USE pegrid
859       USE transpose_indices
860
861       IMPLICIT NONE
862
863       character(len=3) ::  myth_char
864
865       INTEGER ::  i, j, k, m, n, omp_get_thread_num, tn
866
867       REAL, DIMENSION(0:nx)                            ::  work_fftx
868       REAL, DIMENSION(0:nx,1:nz)                       ::  work_trix
869       REAL, DIMENSION(nnx,1:nza,nys_x:nyn_xa,pdims(1)) ::  ar
870       REAL, DIMENSION(:,:,:,:), ALLOCATABLE            ::  tri
871
872
873       CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'start' )
874
875       ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
876
877       tn = 0              ! Default thread number in case of one thread
878!$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
879       DO  j = nys_x, nyn_x
880
881!$        tn = omp_get_thread_num()
882
883          IF ( host(1:3) == 'nec' )  THEN
884!
885!--          Code optimized for vector processors
886             DO  k = 1, nz
887
888                m = 0
889                DO  n = 1, pdims(1)
890                   DO  i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
891                      work_trix(m,k) = ar(i,k,j,n)
892                      m = m + 1
893                   ENDDO
894                ENDDO
895
896             ENDDO
897
898             CALL fft_x_m( work_trix, 'forward' )
899
900          ELSE
901!
902!--          Cache optimized code
903             DO  k = 1, nz
904
905                m = 0
906                DO  n = 1, pdims(1)
907                   DO  i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
908                      work_fftx(m) = ar(i,k,j,n)
909                      m = m + 1
910                   ENDDO
911                ENDDO
912
913                CALL fft_x( work_fftx, 'forward' )
914
915                DO  i = 0, nx
916                   work_trix(i,k) = work_fftx(i)
917                ENDDO
918
919             ENDDO
920
921          ENDIF
922
923!
924!--       Solve the linear equation system
925          CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
926
927          IF ( host(1:3) == 'nec' )  THEN
928!
929!--          Code optimized for vector processors
930             CALL fft_x_m( work_trix, 'backward' )
931
932             DO  k = 1, nz
933
934                m = 0
935                DO  n = 1, pdims(1)
936                   DO  i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
937                      ar(i,k,j,n) = work_trix(m,k)
938                      m = m + 1
939                   ENDDO
940                ENDDO
941
942             ENDDO
943
944          ELSE
945!
946!--          Cache optimized code
947             DO  k = 1, nz
948
949                DO  i = 0, nx
950                   work_fftx(i) = work_trix(i,k)
951                ENDDO
952
953                CALL fft_x( work_fftx, 'backward' )
954
955                m = 0
956                DO  n = 1, pdims(1)
957                   DO  i = 1, nnx_pe( n-1 ) ! WARN: pcoord(i) should be used!!
958                      ar(i,k,j,n) = work_fftx(m)
959                      m = m + 1
960                   ENDDO
961                ENDDO
962
963             ENDDO
964
965          ENDIF
966
967       ENDDO
968
969       DEALLOCATE( tri )
970
971       CALL cpu_log( log_point_s(33), 'fft_x + tridia', 'stop' )
972
973    END SUBROUTINE fftx_tri_fftx
974
975
976    SUBROUTINE fftx_tr_xy( f_in, work, f_out )
977
978!------------------------------------------------------------------------------!
979!  Fourier-transformation along x with subsequent transposition x --> y for
980!  a 1d-decomposition along y
981!
982!  ATTENTION: The NEC-branch of this routine may significantly profit from
983!             further optimizations. So far, performance is much worse than
984!             for routine ffty_tr_yx (more than three times slower).
985!------------------------------------------------------------------------------!
986
987       USE control_parameters
988       USE cpulog
989       USE indices
990       USE interfaces
991       USE pegrid
992       USE transpose_indices
993
994       IMPLICIT NONE
995
996       INTEGER            ::  i, j, k
997
998       REAL, DIMENSION(0:nx,1:nz,nys:nyn)               ::  work_fftx
999       REAL, DIMENSION(1:nza,nys:nyna,0:nxa)            ::  f_in
1000       REAL, DIMENSION(nny,1:nza,nxl_y:nxr_ya,pdims(2)) ::  f_out
1001       REAL, DIMENSION(nys:nyna,1:nza,0:nxa)            ::  work
1002
1003!
1004!--    Carry out the FFT along x, where all data are present due to the
1005!--    1d-decomposition along y. Resort the data in a way that y becomes
1006!--    the first index.
1007       CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
1008
1009       IF ( host(1:3) == 'nec' )  THEN
1010!
1011!--       Code for vector processors
1012!$OMP     PARALLEL PRIVATE ( i, j, k )
1013!$OMP     DO
1014          DO  i = 0, nx
1015
1016             DO  j = nys, nyn
1017                DO  k = 1, nz
1018                   work_fftx(i,k,j) = f_in(k,j,i)
1019                ENDDO
1020             ENDDO
1021
1022          ENDDO
1023
1024!$OMP     DO
1025          DO  j = nys, nyn
1026
1027             CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1028
1029             DO  k = 1, nz
1030                DO  i = 0, nx
1031                   work(j,k,i) = work_fftx(i,k,j)
1032                ENDDO
1033             ENDDO
1034
1035          ENDDO
1036!$OMP     END PARALLEL
1037
1038       ELSE
1039
1040!
1041!--       Cache optimized code (there might be still a potential for better
1042!--       optimization).
1043!$OMP     PARALLEL PRIVATE (i,j,k,work_fftx)
1044!$OMP     DO
1045          DO  i = 0, nx
1046
1047             DO  j = nys, nyn
1048                DO  k = 1, nz
1049                   work_fftx(i,k,j) = f_in(k,j,i)
1050                ENDDO
1051             ENDDO
1052
1053          ENDDO
1054
1055!$OMP     DO
1056          DO  j = nys, nyn
1057             DO  k = 1, nz
1058
1059                CALL fft_x( work_fftx(0:nx,k,j), 'forward' )
1060
1061                DO  i = 0, nx
1062                   work(j,k,i) = work_fftx(i,k,j)
1063                ENDDO
1064             ENDDO
1065
1066          ENDDO
1067!$OMP     END PARALLEL
1068
1069       ENDIF
1070       CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
1071
1072!
1073!--    Transpose array
1074       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1075       CALL MPI_ALLTOALL( work(nys,1,0),      sendrecvcount_xy, MPI_REAL, &
1076                          f_out(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1077                          comm1dy, ierr )
1078       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1079
1080    END SUBROUTINE fftx_tr_xy
1081
1082
1083    SUBROUTINE tr_yx_fftx( f_in, work, f_out )
1084
1085!------------------------------------------------------------------------------!
1086!  Transposition y --> x with a subsequent backward Fourier transformation for
1087!  a 1d-decomposition along x
1088!------------------------------------------------------------------------------!
1089
1090       USE control_parameters
1091       USE cpulog
1092       USE indices
1093       USE interfaces
1094       USE pegrid
1095       USE transpose_indices
1096
1097       IMPLICIT NONE
1098
1099       INTEGER            ::  i, j, k
1100
1101       REAL, DIMENSION(0:nx,1:nz,nys:nyn)               ::  work_fftx
1102       REAL, DIMENSION(nny,1:nza,nxl_y:nxr_ya,pdims(2)) ::  f_in
1103       REAL, DIMENSION(1:nza,nys:nyna,0:nxa)            ::  f_out
1104       REAL, DIMENSION(nys:nyna,1:nza,0:nxa)            ::  work
1105
1106!
1107!--    Transpose array
1108       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1109       CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, &
1110                          work(nys,1,0),     sendrecvcount_xy, MPI_REAL, &
1111                          comm1dy, ierr )
1112       CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1113
1114!
1115!--    Carry out the FFT along x, where all data are present due to the
1116!--    1d-decomposition along y. Resort the data in a way that y becomes
1117!--    the first index.
1118       CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
1119
1120       IF ( host(1:3) == 'nec' )  THEN
1121!
1122!--       Code optimized for vector processors
1123!$OMP     PARALLEL PRIVATE ( i, j, k )
1124!$OMP     DO
1125          DO  j = nys, nyn
1126
1127             DO  k = 1, nz
1128                DO  i = 0, nx
1129                   work_fftx(i,k,j) = work(j,k,i)
1130                ENDDO
1131             ENDDO
1132
1133             CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1134
1135          ENDDO
1136
1137!$OMP     DO
1138          DO  i = 0, nx
1139             DO  j = nys, nyn
1140                DO  k = 1, nz
1141                   f_out(k,j,i) = work_fftx(i,k,j)
1142                ENDDO
1143             ENDDO
1144          ENDDO
1145!$OMP     END PARALLEL
1146
1147       ELSE
1148
1149!
1150!--       Cache optimized code (there might be still a potential for better
1151!--       optimization).
1152!$OMP     PARALLEL PRIVATE (i,j,k,work_fftx)
1153!$OMP     DO
1154          DO  j = nys, nyn
1155             DO  k = 1, nz
1156
1157                DO  i = 0, nx
1158                   work_fftx(i,k,j) = work(j,k,i)
1159                ENDDO
1160
1161                CALL fft_x( work_fftx(0:nx,k,j), 'backward' )
1162
1163             ENDDO
1164          ENDDO
1165
1166!$OMP     DO
1167          DO  i = 0, nx
1168             DO  j = nys, nyn
1169                DO  k = 1, nz
1170                   f_out(k,j,i) = work_fftx(i,k,j)
1171                ENDDO
1172             ENDDO
1173          ENDDO
1174!$OMP     END PARALLEL
1175
1176       ENDIF
1177       CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
1178
1179    END SUBROUTINE tr_yx_fftx
1180
1181
1182    SUBROUTINE ffty_tri_ffty( ar )
1183
1184!------------------------------------------------------------------------------!
1185!  FFT along y, solution of the tridiagonal system and backward FFT for
1186!  a 1d-decomposition along y
1187!
1188!  WARNING: this subroutine may still not work for hybrid parallelization
1189!           with OpenMP (for possible necessary changes see the original
1190!           routine poisfft_hybrid, developed by Klaus Ketelsen, May 2002)
1191!------------------------------------------------------------------------------!
1192
1193       USE control_parameters
1194       USE cpulog
1195       USE grid_variables
1196       USE indices
1197       USE interfaces
1198       USE pegrid
1199       USE transpose_indices
1200
1201       IMPLICIT NONE
1202
1203       INTEGER ::  i, j, k, m, n, omp_get_thread_num, tn
1204
1205       REAL, DIMENSION(0:ny)                            ::  work_ffty
1206       REAL, DIMENSION(0:ny,1:nz)                       ::  work_triy
1207       REAL, DIMENSION(nny,1:nza,nxl_y:nxr_ya,pdims(2)) ::  ar
1208       REAL, DIMENSION(:,:,:,:), ALLOCATABLE            ::  tri
1209
1210
1211       CALL cpu_log( log_point_s(39), 'fft_y + tridia', 'start' )
1212
1213       ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1214
1215       tn = 0           ! Default thread number in case of one thread
1216!$OMP  PARALLEL PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1217!$OMP  DO
1218       DO  i = nxl_y, nxr_y
1219
1220!$        tn = omp_get_thread_num()
1221
1222          IF ( host(1:3) == 'nec' )  THEN
1223!
1224!--          Code optimized for vector processors
1225             DO  k = 1, nz
1226
1227                m = 0
1228                DO  n = 1, pdims(2)
1229                   DO  j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1230                      work_triy(m,k) = ar(j,k,i,n)
1231                      m = m + 1
1232                   ENDDO
1233                ENDDO
1234
1235             ENDDO
1236
1237             CALL fft_y_m( work_triy, ny, 'forward' )
1238
1239          ELSE
1240!
1241!--          Cache optimized code
1242             DO  k = 1, nz
1243
1244                m = 0
1245                DO  n = 1, pdims(2)
1246                   DO  j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1247                      work_ffty(m) = ar(j,k,i,n)
1248                      m = m + 1
1249                   ENDDO
1250                ENDDO
1251
1252                CALL fft_y( work_ffty, 'forward' )
1253
1254                DO  j = 0, ny
1255                   work_triy(j,k) = work_ffty(j)
1256                ENDDO
1257
1258             ENDDO
1259
1260          ENDIF
1261
1262!
1263!--       Solve the linear equation system
1264          CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1265
1266          IF ( host(1:3) == 'nec' )  THEN
1267!
1268!--          Code optimized for vector processors
1269             CALL fft_y_m( work_triy, ny, 'backward' )
1270
1271             DO  k = 1, nz
1272
1273                m = 0
1274                DO  n = 1, pdims(2)
1275                   DO  j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1276                      ar(j,k,i,n) = work_triy(m,k)
1277                      m = m + 1
1278                   ENDDO
1279                ENDDO
1280
1281             ENDDO
1282
1283          ELSE
1284!
1285!--          Cache optimized code
1286             DO  k = 1, nz
1287
1288                DO  j = 0, ny
1289                   work_ffty(j) = work_triy(j,k)
1290                ENDDO
1291
1292                CALL fft_y( work_ffty, 'backward' )
1293
1294                m = 0
1295                DO  n = 1, pdims(2)
1296                   DO  j = 1, nny_pe( n-1 ) ! WARN: pcoord(j) should be used!!
1297                      ar(j,k,i,n) = work_ffty(m)
1298                      m = m + 1
1299                   ENDDO
1300                ENDDO
1301
1302             ENDDO
1303
1304          ENDIF
1305
1306       ENDDO
1307!$OMP  END PARALLEL
1308
1309       DEALLOCATE( tri )
1310
1311       CALL cpu_log( log_point_s(39), 'fft_y + tridia', 'stop' )
1312
1313    END SUBROUTINE ffty_tri_ffty
1314
1315
1316    SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri )
1317
1318!------------------------------------------------------------------------------!
1319! Solves the linear system of equations for a 1d-decomposition along x (see
1320! tridia)
1321!
1322! Attention:  when using the intel compiler, array tri must be passed as an
1323!             argument to the contained subroutines. Otherwise addres faults
1324!             will occur.
1325!             On NEC, tri should not be passed (except for routine substi_1dd)
1326!             because this causes very bad performance.
1327!------------------------------------------------------------------------------!
1328
1329       USE arrays_3d
1330       USE control_parameters
1331
1332       USE pegrid
1333
1334       IMPLICIT NONE
1335
1336       INTEGER ::  i, j, k, nnyh, nx, ny, omp_get_thread_num, tn
1337
1338       REAL    ::  ddx2, ddy2
1339
1340       REAL, DIMENSION(0:nx,1:nz)     ::  ar
1341       REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
1342       REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
1343
1344
1345       nnyh = ( ny + 1 ) / 2
1346
1347!
1348!--    Define constant elements of the tridiagonal matrix.
1349!--    The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
1350!--    the exchanged loops create bank conflicts. The following directive
1351!--    prohibits loop exchange and the loops perform much better.
1352!       tn = omp_get_thread_num()
1353!       WRITE( 120+tn, * ) '+++ id=',myid,' nx=',nx,' thread=', omp_get_thread_num()
1354!       CALL local_flush( 120+tn )
1355!CDIR NOLOOPCHG
1356       DO  k = 0, nz-1
1357          DO  i = 0,nx
1358             tri(2,i,k) = ddzu(k+1) * ddzw(k+1)
1359             tri(3,i,k) = ddzu(k+2) * ddzw(k+1)
1360          ENDDO
1361       ENDDO
1362!       WRITE( 120+tn, * ) '+++ id=',myid,' end of first tridia loop   thread=', omp_get_thread_num()
1363!       CALL local_flush( 120+tn )
1364
1365       IF ( j <= nnyh )  THEN
1366#if defined( __lc )
1367          CALL maketri_1dd( j, tri )
1368#else
1369          CALL maketri_1dd( j )
1370#endif
1371       ELSE
1372#if defined( __lc )
1373          CALL maketri_1dd( ny+1-j, tri )
1374#else
1375          CALL maketri_1dd( ny+1-j )
1376#endif
1377       ENDIF
1378#if defined( __lc )
1379       CALL split_1dd( tri )
1380#else
1381       CALL split_1dd
1382#endif
1383       CALL substi_1dd( ar, tri )
1384
1385    CONTAINS
1386
1387#if defined( __lc )
1388       SUBROUTINE maketri_1dd( j, tri )
1389#else
1390       SUBROUTINE maketri_1dd( j )
1391#endif
1392
1393!------------------------------------------------------------------------------!
1394! computes the i- and j-dependent component of the matrix
1395!------------------------------------------------------------------------------!
1396
1397          USE constants
1398
1399          IMPLICIT NONE
1400
1401          INTEGER ::  i, j, k, nnxh
1402          REAL    ::  a, c
1403
1404          REAL, DIMENSION(0:nx) ::  l
1405
1406#if defined( __lc )
1407          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
1408#endif
1409
1410
1411          nnxh = ( nx + 1 ) / 2
1412!
1413!--       Provide the tridiagonal matrix for solution of the Poisson equation in
1414!--       Fourier space. The coefficients are computed following the method of
1415!--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
1416!--       Siano's original version by discretizing the Poisson equation,
1417!--       before it is Fourier-transformed
1418          DO  i = 0, nx
1419             IF ( i >= 0 .AND. i <= nnxh ) THEN
1420                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * i ) / &
1421                                         FLOAT( nx+1 ) ) ) * ddx2 + &
1422                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
1423                                         FLOAT( ny+1 ) ) ) * ddy2
1424             ELSE
1425                l(i) = 2.0 * ( 1.0 - COS( ( 2.0 * pi * ( nx+1-i ) ) / &
1426                                         FLOAT( nx+1 ) ) ) * ddx2 + &
1427                       2.0 * ( 1.0 - COS( ( 2.0 * pi * j ) / &
1428                                         FLOAT( ny+1 ) ) ) * ddy2
1429             ENDIF
1430          ENDDO
1431
1432          DO  k = 0, nz-1
1433             DO  i = 0, nx
1434                a = -1.0 * ddzu(k+2) * ddzw(k+1)
1435                c = -1.0 * ddzu(k+1) * ddzw(k+1)
1436                tri(1,i,k) = a + c - l(i)
1437             ENDDO
1438          ENDDO
1439          IF ( ibc_p_b == 1  .OR.  ibc_p_b == 2 )  THEN
1440             DO  i = 0, nx
1441                tri(1,i,0) = tri(1,i,0) + tri(2,i,0)
1442             ENDDO
1443          ENDIF
1444          IF ( ibc_p_t == 1 )  THEN
1445             DO  i = 0, nx
1446                tri(1,i,nz-1) = tri(1,i,nz-1) + tri(3,i,nz-1)
1447             ENDDO
1448          ENDIF
1449
1450       END SUBROUTINE maketri_1dd
1451
1452
1453#if defined( __lc )
1454       SUBROUTINE split_1dd( tri )
1455#else
1456       SUBROUTINE split_1dd
1457#endif
1458
1459!------------------------------------------------------------------------------!
1460! Splitting of the tridiagonal matrix (Thomas algorithm)
1461!------------------------------------------------------------------------------!
1462
1463          IMPLICIT NONE
1464
1465          INTEGER ::  i, k
1466
1467#if defined( __lc )
1468          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
1469#endif
1470
1471
1472!
1473!--       Splitting
1474          DO  i = 0, nx
1475             tri(4,i,0) = tri(1,i,0)
1476          ENDDO
1477          DO  k = 1, nz-1
1478             DO  i = 0, nx
1479                tri(5,i,k) = tri(2,i,k) / tri(4,i,k-1)
1480                tri(4,i,k) = tri(1,i,k) - tri(3,i,k-1) * tri(5,i,k)
1481             ENDDO
1482          ENDDO
1483
1484       END SUBROUTINE split_1dd
1485
1486
1487       SUBROUTINE substi_1dd( ar, tri )
1488
1489!------------------------------------------------------------------------------!
1490! Substitution (Forward and Backward) (Thomas algorithm)
1491!------------------------------------------------------------------------------!
1492
1493          IMPLICIT NONE
1494
1495          INTEGER ::  i, k
1496
1497          REAL, DIMENSION(0:nx,nz)       ::  ar
1498          REAL, DIMENSION(0:nx,0:nz-1)   ::  ar1
1499          REAL, DIMENSION(5,0:nx,0:nz-1) ::  tri
1500
1501!
1502!--       Forward substitution
1503          DO  i = 0, nx
1504             ar1(i,0) = ar(i,1)
1505          ENDDO
1506          DO  k = 1, nz-1
1507             DO  i = 0, nx
1508                ar1(i,k) = ar(i,k+1) - tri(5,i,k) * ar1(i,k-1)
1509             ENDDO
1510          ENDDO
1511
1512!
1513!--       Backward substitution
1514          DO  i = 0, nx
1515             ar(i,nz) = ar1(i,nz-1) / tri(4,i,nz-1)
1516          ENDDO
1517          DO  k = nz-2, 0, -1
1518             DO  i = 0, nx
1519                ar(i,k+1) = ( ar1(i,k) - tri(3,i,k) * ar(i,k+2) ) &
1520                            / tri(4,i,k)
1521             ENDDO
1522          ENDDO
1523
1524!
1525!--       Indices i=0, j=0 correspond to horizontally averaged pressure.
1526!--       The respective values of ar should be zero at all k-levels if
1527!--       acceleration of horizontally averaged vertical velocity is zero.
1528          IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
1529             IF ( j == 0 )  THEN
1530                DO  k = 1, nz
1531                   ar(0,k) = 0.0
1532                ENDDO
1533             ENDIF
1534          ENDIF
1535
1536       END SUBROUTINE substi_1dd
1537
1538    END SUBROUTINE tridia_1dd
1539
1540#endif
1541
1542 END MODULE poisfft_mod
Note: See TracBrowser for help on using the repository browser.