source: palm/trunk/SOURCE/poisfft_mod.f90 @ 4680

Last change on this file since 4680 was 4671, checked in by pavelkrc, 4 years ago

Radiative transfer model RTM version 4.1

  • Property svn:keywords set to Id
File size: 45.5 KB
Line 
1!> @file poisfft_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: poisfft_mod.f90 4671 2020-09-09 20:27:58Z gronemeier $
27! OMP bugfix
28!
29! 4649 2020-08-25 12:11:17Z raasch
30! File re-formatted to follow the PALM coding standard
31!
32! 4429 2020-02-27 15:24:30Z raasch
33! Statements added to avoid compile errors due to unused dummy arguments in serial mode
34!
35! 4366 2020-01-09 08:12:43Z raasch
36! Modification concerning NEC vectorizatio
37!
38! 4360 2020-01-07 11:25:50Z suehring
39! Corrected "Former revisions" section
40!
41! 3690 2019-01-22 22:56:42Z knoop
42! OpenACC port for SPEC
43!
44! Revision 1.1  1997/07/24 11:24:14  raasch
45! Initial revision
46!
47!
48!--------------------------------------------------------------------------------------------------!
49! Description:
50! ------------
51!> Solves the Poisson equation with a 2D spectral method
52!>        d^2 p / dx^2 + d^2 p / dy^2 + d^2 p / dz^2 = s
53!>
54!> Input:
55!> real   ar   contains (nnz,nny,nnx) elements of the velocity divergence, starting from (1,nys,nxl)
56!>
57!> Output:
58!> real   ar   contains the solution for perturbation pressure p
59!--------------------------------------------------------------------------------------------------!
60 MODULE poisfft_mod
61
62
63    USE fft_xy,                                                                                    &
64        ONLY:  fft_init,                                                                           &
65               fft_y,                                                                              &
66               fft_y_1d,                                                                           &
67               fft_y_m,                                                                            &
68               fft_x,                                                                              &
69               fft_x_1d,                                                                           &
70               fft_x_m,                                                                            &
71               temperton_fft_vec
72
73    USE indices,                                                                                   &
74        ONLY:  nnx,                                                                                &
75               nny,                                                                                &
76               nx,                                                                                 &
77               nxl,                                                                                &
78               nxr,                                                                                &
79               ny,                                                                                 &
80               nys,                                                                                &
81               nyn,                                                                                &
82               nz
83
84    USE transpose_indices,                                                                         &
85        ONLY:  nxl_y,                                                                              &
86               nxl_z,                                                                              &
87               nxr_y,                                                                              &
88               nxr_z,                                                                              &
89               nys_x,                                                                              &
90               nys_z,                                                                              &
91               nyn_x,                                                                              &
92               nyn_z,                                                                              &
93               nzb_x,                                                                              &
94               nzb_y,                                                                              &
95               nzt_x,                                                                              &
96               nzt_y
97
98    USE tridia_solver,                                                                             &
99        ONLY:  tridia_1dd,                                                                         &
100               tridia_init,                                                                        &
101               tridia_substi,                                                                      &
102               tridia_substi_overlap
103
104    IMPLICIT NONE
105
106    LOGICAL, SAVE ::  poisfft_initialized = .FALSE.  !<
107
108    PRIVATE
109
110    PUBLIC  poisfft, poisfft_init
111
112    INTERFACE poisfft
113       MODULE PROCEDURE poisfft
114    END INTERFACE poisfft
115
116    INTERFACE poisfft_init
117       MODULE PROCEDURE poisfft_init
118    END INTERFACE poisfft_init
119
120
121 CONTAINS
122
123!--------------------------------------------------------------------------------------------------!
124! Description:
125! ------------
126!> Setup coefficients for FFT and the tridiagonal solver
127!--------------------------------------------------------------------------------------------------!
128 SUBROUTINE poisfft_init
129
130    IMPLICIT NONE
131
132
133    CALL fft_init
134
135    CALL tridia_init
136
137    poisfft_initialized = .TRUE.
138
139 END SUBROUTINE poisfft_init
140
141
142
143!--------------------------------------------------------------------------------------------------!
144! Description:
145! ------------
146!> Two-dimensional Fourier Transformation in x- and y-direction.
147!--------------------------------------------------------------------------------------------------!
148 SUBROUTINE poisfft( ar )
149
150    USE control_parameters,                                                                        &
151        ONLY:  transpose_compute_overlap
152
153    USE cpulog,                                                                                    &
154        ONLY:  cpu_log,                                                                            &
155               cpu_log_nowait,                                                                     &
156               log_point_s
157
158    USE kinds
159
160    USE pegrid
161
162    IMPLICIT NONE
163
164    INTEGER(iwp) ::  ii           !<
165    INTEGER(iwp) ::  iind         !<
166    INTEGER(iwp) ::  inew         !<
167    INTEGER(iwp) ::  jj           !<
168    INTEGER(iwp) ::  jind         !<
169    INTEGER(iwp) ::  jnew         !<
170    INTEGER(iwp) ::  ki           !<
171    INTEGER(iwp) ::  kk           !<
172    INTEGER(iwp) ::  knew         !<
173    INTEGER(iwp) ::  n            !<
174    INTEGER(iwp) ::  nblk         !<
175    INTEGER(iwp) ::  nnx_y        !<
176    INTEGER(iwp) ::  nny_z        !<
177    INTEGER(iwp) ::  nnz_x        !<
178    INTEGER(iwp) ::  nxl_y_bound  !<
179    INTEGER(iwp) ::  nxr_y_bound  !<
180
181    INTEGER(iwp), DIMENSION(4) ::  isave  !<
182
183    REAL(wp), DIMENSION(1:nz,nys:nyn,nxl:nxr) ::  ar      !<
184    REAL(wp), DIMENSION(nys:nyn,nxl:nxr,1:nz) ::  ar_inv  !<
185
186#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
187#if __acc_fft_device
188    !$ACC DECLARE CREATE(ar_inv)
189#endif
190
191    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  ar1      !<
192    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_in     !<
193    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_inv    !<
194    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_out_y  !<
195    REAL(wp), DIMENSION(:,:,:),   ALLOCATABLE ::  f_out_z  !<
196
197
198    CALL cpu_log( log_point_s(3), 'poisfft', 'start' )
199
200    IF ( .NOT. poisfft_initialized )  CALL poisfft_init
201
202#if !__acc_fft_device
203    !$ACC UPDATE HOST(ar)
204#endif
205
206#ifndef _OPENACC
207!
208!-- Two-dimensional Fourier Transformation in x- and y-direction.
209    IF ( pdims(2) == 1  .AND.  pdims(1) > 1 )  THEN
210
211!
212!--    1d-domain-decomposition along x:
213!--    FFT along y and transposition y --> x
214       CALL ffty_tr_yx( ar, ar )
215
216!
217!--    FFT along x, solving the tridiagonal system and backward FFT
218       CALL fftx_tri_fftx( ar )
219
220!
221!--    Transposition x --> y and backward FFT along y
222       CALL tr_xy_ffty( ar, ar )
223
224    ELSEIF ( pdims(1) == 1 .AND. pdims(2) > 1 )  THEN
225
226!
227!--    1d-domain-decomposition along y:
228!--    FFT along x and transposition x --> y
229       CALL fftx_tr_xy( ar, ar )
230
231!
232!--    FFT along y, solving the tridiagonal system and backward FFT
233       CALL ffty_tri_ffty( ar )
234
235!
236!--    Transposition y --> x and backward FFT along x
237       CALL tr_yx_fftx( ar, ar )
238
239    ELSEIF ( .NOT. transpose_compute_overlap )  THEN
240#endif
241
242!
243!--    2d-domain-decomposition or no decomposition (1 PE run)
244!--    Transposition z --> x
245       CALL cpu_log( log_point_s(5), 'transpo forward', 'start' )
246       CALL resort_for_zx( ar, ar_inv )
247       CALL transpose_zx( ar_inv, ar )
248       CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
249
250       CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
251       IF ( temperton_fft_vec )  THEN
252!
253!--       Vector version outputs a transformed array ar_inv that does not require resorting
254!--       (which is done for ar further below)
255          CALL fft_x( ar, 'forward',  ar_inv=ar_inv)
256       ELSE
257          CALL fft_x( ar, 'forward')
258       ENDIF
259       CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
260
261!
262!--    Transposition x --> y
263       CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
264       IF( .NOT. temperton_fft_vec )  CALL resort_for_xy( ar, ar_inv )
265       CALL transpose_xy( ar_inv, ar )
266       CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
267
268       CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
269       IF ( temperton_fft_vec )  THEN
270!
271!--       Input array ar_inv from fft_x can be directly used here.
272!--       The output (also in array ar_inv) does not require resorting below.
273          CALL fft_y( ar, 'forward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,    &
274                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
275       ELSE
276          CALL fft_y( ar, 'forward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,         &
277                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
278       ENDIF
279       CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
280
281!
282!--    Transposition y --> z
283       CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
284       IF ( .NOT. temperton_fft_vec )  CALL resort_for_yz( ar, ar_inv )
285       CALL transpose_yz( ar_inv, ar )
286       CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
287
288!
289!--    Solve the tridiagonal equation system along z
290       CALL cpu_log( log_point_s(6), 'tridia', 'start' )
291       CALL tridia_substi( ar )
292       CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
293
294!
295!--    Inverse Fourier Transformation
296!--    Transposition z --> y
297       CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
298       CALL transpose_zy( ar, ar_inv )
299!
300!--    The fft_y below (vector branch) can directly process ar_inv (i.e. does not require a
301!--    resorting)
302       IF ( .NOT. temperton_fft_vec )  CALL resort_for_zy( ar_inv, ar )
303       CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
304
305       CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
306       IF ( temperton_fft_vec )  THEN
307!
308!--       Output array ar_inv can be used as input to the below fft_x routine without resorting
309          CALL fft_y( ar, 'backward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,   &
310                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
311       ELSE
312          CALL fft_y( ar, 'backward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,        &
313                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
314       ENDIF
315
316       CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
317
318!
319!--    Transposition y --> x
320       CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
321       CALL transpose_yx( ar, ar_inv )
322       IF ( .NOT. temperton_fft_vec )  CALL resort_for_yx( ar_inv, ar )
323       CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
324
325       CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
326       IF ( temperton_fft_vec )  THEN
327          CALL fft_x( ar, 'backward',  ar_inv=ar_inv )
328       ELSE
329          CALL fft_x( ar, 'backward' )
330       ENDIF
331       CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
332
333!
334!--    Transposition x --> z
335       CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
336       CALL transpose_xz( ar, ar_inv )
337       CALL resort_for_xz( ar_inv, ar )
338       CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
339
340#ifndef _OPENACC
341    ELSE
342
343!
344!--    2d-domain-decomposition or no decomposition (1 PE run) with overlapping transposition / fft
345!--    cputime logging must not use barriers, which would prevent overlapping
346       ALLOCATE( f_out_y(0:ny,nxl_y:nxr_y,nzb_y:nzt_y), &
347                 f_out_z(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
348!
349!--    Transposition z --> x + subsequent fft along x
350       ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
351       CALL resort_for_zx( ar, f_inv )
352!
353!--    Save original indices and gridpoint counter
354       isave(1) = nz
355       isave(2) = nzb_x
356       isave(3) = nzt_x
357       isave(4) = sendrecvcount_zx
358!
359!--    Set new indices for transformation
360       nblk  = nz / pdims(1)
361       nz    = pdims(1)
362       nnz_x = 1
363       nzb_x = 1 + myidx * nnz_x
364       nzt_x = ( myidx + 1 ) * nnz_x
365       sendrecvcount_zx = nnx * nny * nnz_x
366
367       ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
368       ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
369
370       DO  kk = 1, nblk
371
372          IF ( kk == 1 )  THEN
373             CALL cpu_log( log_point_s(5), 'transpo forward', 'start', cpu_log_nowait )
374          ELSE
375             CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
376          ENDIF
377
378          DO  knew = 1, nz
379             ki = kk + nblk * ( knew - 1 )
380             f_in(:,:,knew) = f_inv(:,:,ki)
381          ENDDO
382
383          CALL transpose_zx( f_in, ar1(:,:,:))
384          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
385
386          IF ( kk == 1 )  THEN
387             CALL cpu_log( log_point_s(4), 'fft_x', 'start', cpu_log_nowait )
388          ELSE
389             CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
390          ENDIF
391
392          n = isave(2) + kk - 1
393          CALL fft_x( ar1(:,:,:), 'forward',  ar_2d = f_out_z(:,:,n) )
394          CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
395
396       ENDDO
397!
398!--    Restore original indices/counters
399       nz               = isave(1)
400       nzb_x            = isave(2)
401       nzt_x            = isave(3)
402       sendrecvcount_zx = isave(4)
403
404       DEALLOCATE( ar1, f_in, f_inv )
405
406!
407!--    Transposition x --> y + subsequent fft along y
408       ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
409       CALL resort_for_xy( f_out_z, f_inv )
410!
411!--    Save original indices and gridpoint counter
412       isave(1) = nx
413       isave(2) = nxl_y
414       isave(3) = nxr_y
415       isave(4) = sendrecvcount_xy
416!
417!--    Set new indices for transformation
418       nblk  = ( ( nx+1 ) / pdims(2) ) - 1
419       nx    = pdims(2)
420       nnx_y = 1
421       nxl_y = myidy * nnx_y
422       nxr_y = ( myidy + 1 ) * nnx_y - 1
423       sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
424
425       ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
426       ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
427
428       DO  ii = 0, nblk
429
430          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
431
432          DO  inew = 0, nx-1
433             iind = ii + ( nblk + 1 ) * inew
434             f_in(:,:,inew) = f_inv(:,:,iind)
435          ENDDO
436
437          CALL transpose_xy( f_in, ar1(:,:,:) )
438
439          CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
440
441          IF ( ii == 1 )  THEN
442             CALL cpu_log( log_point_s(7), 'fft_y', 'start', cpu_log_nowait )
443          ELSE
444             CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
445          ENDIF
446
447          nxl_y_bound = isave(2)
448          nxr_y_bound = isave(3)
449          n           = isave(2) + ii
450          CALL fft_y( ar1(:,:,:), 'forward', ar_tr = f_out_y, nxl_y_bound = nxl_y_bound,           &
451                      nxr_y_bound = nxr_y_bound, nxl_y_l = n, nxr_y_l = n )
452
453          CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
454
455       ENDDO
456!
457!--    Restore original indices/counters
458       nx               = isave(1)
459       nxl_y            = isave(2)
460       nxr_y            = isave(3)
461       sendrecvcount_xy = isave(4)
462
463       DEALLOCATE( ar1, f_in, f_inv )
464
465!
466!--    Transposition y --> z + subsequent tridia + resort for z --> y
467       ALLOCATE( f_inv(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
468       CALL resort_for_yz( f_out_y, f_inv )
469!
470!--    Save original indices and gridpoint counter
471       isave(1) = ny
472       isave(2) = nys_z
473       isave(3) = nyn_z
474       isave(4) = sendrecvcount_yz
475!
476!--    Set new indices for transformation
477       nblk             = ( ( ny+1 ) / pdims(1) ) - 1
478       ny               = pdims(1)
479       nny_z            = 1
480       nys_z            = myidx * nny_z
481       nyn_z            = ( myidx + 1 ) * nny_z - 1
482       sendrecvcount_yz = ( nxr_y-nxl_y+1 ) * nny_z * ( nzt_y-nzb_y+1 )
483
484       ALLOCATE( ar1(nxl_z:nxr_z,nys_z:nyn_z,1:nz) )
485       ALLOCATE( f_in(nxl_y:nxr_y,nzb_y:nzt_y,0:ny) )
486
487       DO  jj = 0, nblk
488!
489!--       Forward Fourier Transformation
490!--       Transposition y --> z
491          CALL cpu_log( log_point_s(5), 'transpo forward', 'continue', cpu_log_nowait )
492
493          DO  jnew = 0, ny-1
494             jind = jj + ( nblk + 1 ) * jnew
495             f_in(:,:,jnew) = f_inv(:,:,jind)
496          ENDDO
497
498          CALL transpose_yz( f_in, ar1(:,:,:) )
499
500          IF ( jj == nblk )  THEN
501             CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
502          ELSE
503             CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
504          ENDIF
505
506!
507!--       Solve the tridiagonal equation system along z
508          CALL cpu_log( log_point_s(6), 'tridia', 'start', cpu_log_nowait )
509
510          n = isave(2) + jj
511          CALL tridia_substi_overlap( ar1(:,:,:), n )
512
513          CALL cpu_log( log_point_s(6), 'tridia', 'stop' )
514
515!
516!--       Inverse Fourier Transformation
517!--       Transposition z --> y
518!--       Only one thread should call MPI routines, therefore forward and backward tranpose are in
519!--       the same section
520          IF ( jj == 0 )  THEN
521             CALL cpu_log( log_point_s(8), 'transpo invers', 'start', cpu_log_nowait )
522          ELSE
523             CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
524          ENDIF
525
526          CALL transpose_zy( ar1(:,:,:), f_in )
527
528          DO  jnew = 0, ny-1
529             jind = jj + ( nblk + 1 ) * jnew
530             f_inv(:,:,jind) = f_in(:,:,jnew)
531          ENDDO
532
533          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
534
535       ENDDO
536!
537!--    Restore original indices/counters
538       ny               = isave(1)
539       nys_z            = isave(2)
540       nyn_z            = isave(3)
541       sendrecvcount_yz = isave(4)
542
543       CALL resort_for_zy( f_inv, f_out_y )
544
545       DEALLOCATE( ar1, f_in, f_inv )
546
547!
548!--    fft along y backward + subsequent transposition y --> x
549       ALLOCATE( f_inv(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
550!
551!--    Save original indices and gridpoint counter
552       isave(1) = nx
553       isave(2) = nxl_y
554       isave(3) = nxr_y
555       isave(4) = sendrecvcount_xy
556!
557!--    Set new indices for transformation
558       nblk             = ( ( nx+1 ) / pdims(2) ) - 1
559       nx               = pdims(2)
560       nnx_y            = 1
561       nxl_y            = myidy * nnx_y
562       nxr_y            = ( myidy + 1 ) * nnx_y - 1
563       sendrecvcount_xy = nnx_y * ( nyn_x-nys_x+1 ) * ( nzt_x-nzb_x+1 )
564
565       ALLOCATE( ar1(0:ny,nxl_y:nxr_y,nzb_y:nzt_y) )
566       ALLOCATE( f_in(nys_x:nyn_x,nzb_x:nzt_x,0:nx) )
567
568       DO  ii = 0, nblk
569
570          CALL cpu_log( log_point_s(7), 'fft_y', 'continue', cpu_log_nowait )
571
572          n = isave(2) + ii
573          nxl_y_bound = isave(2)
574          nxr_y_bound = isave(3)
575
576          CALL fft_y( ar1(:,:,:), 'backward', ar_tr = f_out_y, nxl_y_bound = nxl_y_bound,          &
577                      nxr_y_bound = nxr_y_bound, nxl_y_l = n, nxr_y_l = n )
578
579          IF ( ii == nblk )  THEN
580             CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
581          ELSE
582             CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
583          ENDIF
584
585          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
586
587          CALL transpose_yx( ar1(:,:,:), f_in )
588
589          DO  inew = 0, nx-1
590             iind = ii + (nblk+1) * inew
591             f_inv(:,:,iind) = f_in(:,:,inew)
592          ENDDO
593
594          CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
595
596       ENDDO
597!
598!--    Restore original indices/counters
599       nx               = isave(1)
600       nxl_y            = isave(2)
601       nxr_y            = isave(3)
602       sendrecvcount_xy = isave(4)
603
604       CALL resort_for_yx( f_inv, f_out_z )
605
606       DEALLOCATE( ar1, f_in, f_inv )
607
608!
609!--    fft along x backward + subsequent final transposition x --> z
610       ALLOCATE( f_inv(nys:nyn,nxl:nxr,1:nz) )
611!
612!--    Save original indices and gridpoint counter
613       isave(1) = nz
614       isave(2) = nzb_x
615       isave(3) = nzt_x
616       isave(4) = sendrecvcount_zx
617!
618!--    Set new indices for transformation
619       nblk             = nz / pdims(1)
620       nz               = pdims(1)
621       nnz_x            = 1
622       nzb_x            = 1 + myidx * nnz_x
623       nzt_x            = ( myidx + 1 ) * nnz_x
624       sendrecvcount_zx = nnx * nny * nnz_x
625
626       ALLOCATE( ar1(0:nx,nys_x:nyn_x,nzb_x:nzt_x) )
627       ALLOCATE( f_in(nys:nyn,nxl:nxr,1:nz) )
628
629       DO  kk = 1, nblk
630
631          CALL cpu_log( log_point_s(4), 'fft_x', 'continue', cpu_log_nowait )
632
633          n = isave(2) + kk - 1
634          CALL fft_x( ar1(:,:,:), 'backward', f_out_z(:,:,n) )
635
636          IF ( kk == nblk )  THEN
637             CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
638          ELSE
639             CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
640          ENDIF
641
642          CALL cpu_log( log_point_s(8), 'transpo invers', 'continue', cpu_log_nowait )
643
644          CALL transpose_xz( ar1(:,:,:), f_in )
645
646          DO  knew = 1, nz
647             ki = kk + nblk * ( knew - 1 )
648             f_inv(:,:,ki) = f_in(:,:,knew)
649          ENDDO
650
651          IF ( kk == nblk )  THEN
652             CALL cpu_log( log_point_s(8), 'transpo invers', 'stop' )
653          ELSE
654             CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
655          ENDIF
656
657       ENDDO
658!
659!--    Restore original indices/counters
660       nz               = isave(1)
661       nzb_x            = isave(2)
662       nzt_x            = isave(3)
663       sendrecvcount_zx = isave(4)
664
665       CALL resort_for_xz( f_inv, ar )
666
667       DEALLOCATE( ar1, f_in, f_inv )
668
669    ENDIF
670#endif
671
672#if !__acc_fft_device
673    !$ACC UPDATE DEVICE(ar)
674#endif
675
676    CALL cpu_log( log_point_s(3), 'poisfft', 'stop' )
677
678 END SUBROUTINE poisfft
679
680
681!--------------------------------------------------------------------------------------------------!
682! Description:
683! ------------
684!> Fourier-transformation along y with subsequent transposition y --> x for a 1d-decomposition
685!> along x.
686!>
687!> @attention The performance of this routine is much faster on the NEC-SX6, if the first index of
688!>            work_ffty_vec is odd. Otherwise memory bank conflicts may occur (especially if the
689!>            index is a multiple of 128). That's why work_ffty_vec is dimensioned as 0:ny+1.
690!>            Of course, this will not work if users are using an odd number of gridpoints along y.
691!--------------------------------------------------------------------------------------------------!
692 SUBROUTINE ffty_tr_yx( f_in, f_out )
693
694    USE control_parameters,                                                                        &
695        ONLY:  loop_optimization
696
697    USE cpulog,                                                                                    &
698        ONLY:  cpu_log,                                                                            &
699               log_point_s
700
701    USE kinds
702
703    USE pegrid
704
705    IMPLICIT NONE
706
707    INTEGER(iwp) ::  i       !<
708    INTEGER(iwp) ::  iend    !<
709    INTEGER(iwp) ::  iouter  !<
710    INTEGER(iwp) ::  ir      !<
711    INTEGER(iwp) ::  j       !<
712    INTEGER(iwp) ::  k       !<
713
714    INTEGER(iwp), PARAMETER ::  stridex = 4  !<
715
716    REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr)             ::  f_in   !<
717    REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  f_out  !<
718    REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny)             ::  work   !<
719
720    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  work_ffty      !<
721    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  work_ffty_vec  !<
722
723!
724!-- Carry out the FFT along y, where all data are present due to the 1d-decomposition along x.
725!-- Resort the data in a way that x becomes the first index.
726    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'start' )
727
728    IF ( loop_optimization == 'vector' )  THEN
729
730       ALLOCATE( work_ffty_vec(0:ny+1,1:nz,nxl:nxr) )
731!
732!--    Code optimized for vector processors
733       !$OMP PARALLEL PRIVATE ( i, j, k )
734       !$OMP DO
735       DO  i = nxl, nxr
736
737          DO  j = 0, ny
738             DO  k = 1, nz
739                work_ffty_vec(j,k,i) = f_in(k,j,i)
740             ENDDO
741          ENDDO
742
743          CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'forward' )
744
745       ENDDO
746
747       !$OMP DO
748       DO  k = 1, nz
749          DO  j = 0, ny
750             DO  i = nxl, nxr
751                work(i,k,j) = work_ffty_vec(j,k,i)
752             ENDDO
753          ENDDO
754       ENDDO
755       !$OMP END PARALLEL
756
757       DEALLOCATE( work_ffty_vec )
758
759    ELSE
760!
761!--    Cache optimized code.
762       ALLOCATE( work_ffty(0:ny,stridex) )
763!
764!--    The i-(x-)direction is split into a strided outer loop and an inner loop for better cache
765!--    performance
766       !$OMP PARALLEL PRIVATE (i,iend,iouter,ir,j,k,work_ffty)
767       !$OMP DO
768       DO  iouter = nxl, nxr, stridex
769
770          iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
771
772          DO  k = 1, nz
773
774             DO  i = iouter, iend
775
776                ir = i-iouter+1  ! Counter within a stride
777                DO  j = 0, ny
778                   work_ffty(j,ir) = f_in(k,j,i)
779                ENDDO
780!
781!--             FFT along y
782                CALL fft_y_1d( work_ffty(:,ir), 'forward' )
783
784             ENDDO
785
786!
787!--          Resort
788             DO  j = 0, ny
789                DO  i = iouter, iend
790                   work(i,k,j) = work_ffty(j,i-iouter+1)
791                ENDDO
792             ENDDO
793
794          ENDDO
795
796       ENDDO
797       !$OMP END PARALLEL
798
799       DEALLOCATE( work_ffty )
800
801    ENDIF
802
803    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'pause' )
804
805!
806!-- Transpose array
807#if defined( __parallel )
808    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
809    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
810    CALL MPI_ALLTOALL( work(nxl,1,0), sendrecvcount_xy, MPI_REAL, f_out(1,1,nys_x,1),              &
811                       sendrecvcount_xy, MPI_REAL, comm1dx, ierr )
812    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
813#else
814!
815!-- Next line required to avoid compile error about unused dummy argument in serial mode
816    i = SIZE( f_out )
817#endif
818
819 END SUBROUTINE ffty_tr_yx
820
821
822!--------------------------------------------------------------------------------------------------!
823! Description:
824! ------------
825!> Transposition x --> y with a subsequent backward Fourier transformation for a 1d-decomposition
826!> along x
827!--------------------------------------------------------------------------------------------------!
828 SUBROUTINE tr_xy_ffty( f_in, f_out )
829
830    USE control_parameters,                                                                        &
831        ONLY:  loop_optimization
832
833    USE cpulog,                                                                                    &
834        ONLY:  cpu_log,                                                                            &
835               log_point_s
836
837    USE kinds
838
839    USE pegrid
840
841    IMPLICIT NONE
842
843    INTEGER(iwp) ::  i       !<
844    INTEGER(iwp) ::  iend    !<
845    INTEGER(iwp) ::  iouter  !<
846    INTEGER(iwp) ::  ir      !<
847    INTEGER(iwp) ::  j       !<
848    INTEGER(iwp) ::  k       !<
849
850    INTEGER(iwp), PARAMETER ::  stridex = 4  !<
851
852    REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  f_in   !<
853    REAL(wp), DIMENSION(1:nz,0:ny,nxl:nxr)             ::  f_out  !<
854    REAL(wp), DIMENSION(nxl:nxr,1:nz,0:ny)             ::  work   !<
855
856    REAL(wp), DIMENSION(:,:),   ALLOCATABLE ::  work_ffty      !<
857    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  work_ffty_vec  !<
858
859!
860!-- Transpose array
861#if defined( __parallel )
862    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
863    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
864    CALL MPI_ALLTOALL( f_in(1,1,nys_x,1), sendrecvcount_xy, MPI_REAL, work(nxl,1,0),               &
865                       sendrecvcount_xy, MPI_REAL, comm1dx, ierr )
866    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
867#else
868!
869!-- Next line required to avoid compile error about unused dummy argument in serial mode
870    i = SIZE( f_in )
871#endif
872
873!
874!-- Resort the data in a way that y becomes the first index and carry out the backward fft along y.
875    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'continue' )
876
877    IF ( loop_optimization == 'vector' )  THEN
878
879       ALLOCATE( work_ffty_vec(0:ny+1,1:nz,nxl:nxr) )
880!
881!--    Code optimized for vector processors
882       !$OMP PARALLEL PRIVATE ( i, j, k )
883       !$OMP DO
884       DO  k = 1, nz
885          DO  j = 0, ny
886             DO  i = nxl, nxr
887                work_ffty_vec(j,k,i) = work(i,k,j)
888             ENDDO
889          ENDDO
890       ENDDO
891
892       !$OMP DO
893       DO  i = nxl, nxr
894
895          CALL fft_y_m( work_ffty_vec(:,:,i), ny+1, 'backward' )
896
897          DO  j = 0, ny
898             DO  k = 1, nz
899                f_out(k,j,i) = work_ffty_vec(j,k,i)
900             ENDDO
901          ENDDO
902
903       ENDDO
904       !$OMP END PARALLEL
905
906       DEALLOCATE( work_ffty_vec )
907
908    ELSE
909!
910!--    Cache optimized code.
911       ALLOCATE( work_ffty(0:ny,stridex) )
912!
913!--    The i-(x-)direction is split into a strided outer loop and an inner loop for better cache
914!--    performance
915       !$OMP PARALLEL PRIVATE ( i, iend, iouter, ir, j, k, work_ffty )
916       !$OMP DO
917       DO  iouter = nxl, nxr, stridex
918
919          iend = MIN( iouter+stridex-1, nxr )  ! Upper bound for inner i loop
920
921          DO  k = 1, nz
922!
923!--          Resort
924             DO  j = 0, ny
925                DO  i = iouter, iend
926                   work_ffty(j,i-iouter+1) = work(i,k,j)
927                ENDDO
928             ENDDO
929
930             DO  i = iouter, iend
931
932!
933!--             FFT along y
934                ir = i-iouter+1  ! Counter within a stride
935                CALL fft_y_1d( work_ffty(:,ir), 'backward' )
936
937                DO  j = 0, ny
938                   f_out(k,j,i) = work_ffty(j,ir)
939                ENDDO
940             ENDDO
941
942          ENDDO
943
944       ENDDO
945       !$OMP END PARALLEL
946
947       DEALLOCATE( work_ffty )
948
949    ENDIF
950
951    CALL cpu_log( log_point_s(7), 'fft_y_1d', 'stop' )
952
953 END SUBROUTINE tr_xy_ffty
954
955
956!--------------------------------------------------------------------------------------------------!
957! Description:
958! ------------
959!> FFT along x, solution of the tridiagonal system and backward FFT for a 1d-decomposition along x
960!>
961!> @warning This subroutine may still not work for hybrid parallelization with OpenMP (for possible
962!>          necessary changes see the original routine poisfft_hybrid, developed by Klaus Ketelsen,
963!>          May 2002)
964!--------------------------------------------------------------------------------------------------!
965 SUBROUTINE fftx_tri_fftx( ar )
966
967    USE control_parameters,                                                                        &
968        ONLY:  loop_optimization
969
970    USE cpulog,                                                                                    &
971        ONLY:  cpu_log,                                                                            &
972               log_point_s
973
974    USE grid_variables,                                                                            &
975        ONLY:  ddx2,                                                                               &
976               ddy2
977
978    USE kinds
979
980    USE pegrid
981
982    IMPLICIT NONE
983
984    INTEGER(iwp) ::  i                   !<
985    INTEGER(iwp) ::  j                   !<
986    INTEGER(iwp) ::  k                   !<
987    INTEGER(iwp) ::  m                   !<
988    INTEGER(iwp) ::  n                   !<
989!$  INTEGER(iwp) ::  omp_get_thread_num  !<
990    INTEGER(iwp) ::  tn                  !<
991
992    REAL(wp), DIMENSION(0:nx)                          ::  work_fftx  !<
993    REAL(wp), DIMENSION(0:nx,1:nz)                     ::  work_trix  !<
994    REAL(wp), DIMENSION(nnx,1:nz,nys_x:nyn_x,pdims(1)) ::  ar         !<
995    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          ::  tri        !<
996
997
998    CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'start' )
999
1000    ALLOCATE( tri(5,0:nx,0:nz-1,0:threads_per_task-1) )
1001
1002    tn = 0              ! Default thread number in case of one thread
1003!$OMP  PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_fftx, work_trix )
1004    DO  j = nys_x, nyn_x
1005
1006!$     tn = omp_get_thread_num()
1007
1008       IF ( loop_optimization == 'vector' )  THEN
1009!
1010!--       Code optimized for vector processors
1011          DO  k = 1, nz
1012
1013             m = 0
1014             DO  n = 1, pdims(1)
1015                DO  i = 1, nnx
1016                   work_trix(m,k) = ar(i,k,j,n)
1017                   m = m + 1
1018                ENDDO
1019             ENDDO
1020
1021          ENDDO
1022
1023          CALL fft_x_m( work_trix, 'forward' )
1024
1025       ELSE
1026!
1027!--       Cache optimized code
1028          DO  k = 1, nz
1029
1030             m = 0
1031             DO  n = 1, pdims(1)
1032                DO  i = 1, nnx
1033                   work_fftx(m) = ar(i,k,j,n)
1034                   m = m + 1
1035                ENDDO
1036             ENDDO
1037
1038             CALL fft_x_1d( work_fftx, 'forward' )
1039
1040             DO  i = 0, nx
1041                work_trix(i,k) = work_fftx(i)
1042             ENDDO
1043
1044          ENDDO
1045
1046       ENDIF
1047
1048!
1049!--    Solve the linear equation system
1050       CALL tridia_1dd( ddx2, ddy2, nx, ny, j, work_trix, tri(:,:,:,tn) )
1051
1052       IF ( loop_optimization == 'vector' )  THEN
1053!
1054!--       Code optimized for vector processors
1055          CALL fft_x_m( work_trix, 'backward' )
1056
1057          DO  k = 1, nz
1058
1059             m = 0
1060             DO  n = 1, pdims(1)
1061                DO  i = 1, nnx
1062                   ar(i,k,j,n) = work_trix(m,k)
1063                   m = m + 1
1064                ENDDO
1065             ENDDO
1066
1067          ENDDO
1068
1069       ELSE
1070!
1071!--       Cache optimized code
1072          DO  k = 1, nz
1073
1074             DO  i = 0, nx
1075                work_fftx(i) = work_trix(i,k)
1076             ENDDO
1077
1078             CALL fft_x_1d( work_fftx, 'backward' )
1079
1080             m = 0
1081             DO  n = 1, pdims(1)
1082                DO  i = 1, nnx
1083                   ar(i,k,j,n) = work_fftx(m)
1084                   m = m + 1
1085                ENDDO
1086             ENDDO
1087
1088          ENDDO
1089
1090       ENDIF
1091
1092    ENDDO
1093
1094    DEALLOCATE( tri )
1095
1096    CALL cpu_log( log_point_s(33), 'fft_x_1d + tridia', 'stop' )
1097
1098 END SUBROUTINE fftx_tri_fftx
1099
1100
1101!--------------------------------------------------------------------------------------------------!
1102! Description:
1103! ------------
1104!> Fourier-transformation along x with subsequent transposition x --> y for a 1d-decomposition
1105!> along y.
1106!>
1107!> @attention NEC-branch of this routine may significantly profit from further optimizations. So
1108!>            far, performance is much worse than for routine ffty_tr_yx (more than three times
1109!>            slower).
1110!--------------------------------------------------------------------------------------------------!
1111 SUBROUTINE fftx_tr_xy( f_in, f_out )
1112
1113
1114    USE control_parameters,                                                                        &
1115        ONLY:  loop_optimization
1116
1117    USE cpulog,                                                                                    &
1118        ONLY:  cpu_log,                                                                            &
1119               log_point_s
1120
1121    USE kinds
1122
1123    USE pegrid
1124
1125    IMPLICIT NONE
1126
1127    INTEGER(iwp) ::  i  !<
1128    INTEGER(iwp) ::  j  !<
1129    INTEGER(iwp) ::  k  !<
1130
1131    REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn)             ::  work_fftx  !<
1132    REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx)             ::  f_in       !<
1133    REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  f_out      !<
1134    REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx)             ::  work       !<
1135
1136!
1137!--  Carry out the FFT along x, where all data are present due to the 1d-decomposition along y.
1138!--  Resort the data in a way that y becomes the first index.
1139     CALL cpu_log( log_point_s(4), 'fft_x_1d', 'start' )
1140
1141     IF ( loop_optimization == 'vector' )  THEN
1142!
1143!--    Code for vector processors
1144!$OMP  PARALLEL PRIVATE ( i, j, k )
1145!$OMP  DO
1146       DO  i = 0, nx
1147
1148          DO  j = nys, nyn
1149             DO  k = 1, nz
1150                work_fftx(i,k,j) = f_in(k,j,i)
1151             ENDDO
1152          ENDDO
1153
1154       ENDDO
1155
1156!$OMP  DO
1157       DO  j = nys, nyn
1158
1159          CALL fft_x_m( work_fftx(:,:,j), 'forward' )
1160
1161          DO  k = 1, nz
1162             DO  i = 0, nx
1163                work(j,k,i) = work_fftx(i,k,j)
1164             ENDDO
1165          ENDDO
1166
1167       ENDDO
1168!$OMP  END PARALLEL
1169
1170    ELSE
1171
1172!
1173!--    Cache optimized code (there might still be a potential for better optimization).
1174!$OMP  PARALLEL PRIVATE (i,j,k)
1175!$OMP  DO
1176       DO  i = 0, nx
1177
1178          DO  j = nys, nyn
1179             DO  k = 1, nz
1180                work_fftx(i,k,j) = f_in(k,j,i)
1181             ENDDO
1182          ENDDO
1183
1184       ENDDO
1185
1186!$OMP  DO
1187       DO  j = nys, nyn
1188          DO  k = 1, nz
1189
1190             CALL fft_x_1d( work_fftx(0:nx,k,j), 'forward' )
1191
1192             DO  i = 0, nx
1193                work(j,k,i) = work_fftx(i,k,j)
1194             ENDDO
1195          ENDDO
1196
1197       ENDDO
1198!$OMP  END PARALLEL
1199
1200    ENDIF
1201    CALL cpu_log( log_point_s(4), 'fft_x_1d', 'pause' )
1202
1203!
1204!-- Transpose array
1205#if defined( __parallel )
1206    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1207    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1208    CALL MPI_ALLTOALL( work(nys,1,0), sendrecvcount_xy, MPI_REAL, f_out(1,1,nxl_y,1),              &
1209                       sendrecvcount_xy, MPI_REAL, comm1dy, ierr )
1210    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1211#else
1212!
1213!-- Next line required to avoid compile error about unused dummy argument in serial mode
1214    i = SIZE( f_out )
1215#endif
1216
1217 END SUBROUTINE fftx_tr_xy
1218
1219
1220!--------------------------------------------------------------------------------------------------!
1221! Description:
1222! ------------
1223!> Transposition y --> x with a subsequent backward Fourier transformation for a 1d-decomposition
1224!> along x.
1225!--------------------------------------------------------------------------------------------------!
1226 SUBROUTINE tr_yx_fftx( f_in, f_out )
1227
1228
1229    USE control_parameters,                                                                        &
1230        ONLY:  loop_optimization
1231
1232    USE cpulog,                                                                                    &
1233        ONLY:  cpu_log,                                                                            &
1234               log_point_s
1235
1236    USE kinds
1237
1238    USE pegrid
1239
1240    IMPLICIT NONE
1241
1242    INTEGER(iwp) ::  i  !<
1243    INTEGER(iwp) ::  j  !<
1244    INTEGER(iwp) ::  k  !<
1245
1246    REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  f_in       !<
1247    REAL(wp), DIMENSION(1:nz,nys:nyn,0:nx)             ::  f_out      !<
1248    REAL(wp), DIMENSION(nys:nyn,1:nz,0:nx)             ::  work       !<
1249    REAL(wp), DIMENSION(0:nx,1:nz,nys:nyn)             ::  work_fftx  !<
1250
1251
1252!
1253!-- Transpose array
1254#if defined( __parallel )
1255    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'start' )
1256    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1257    CALL MPI_ALLTOALL( f_in(1,1,nxl_y,1), sendrecvcount_xy, MPI_REAL, work(nys,1,0),               &
1258                       sendrecvcount_xy, MPI_REAL, comm1dy, ierr )
1259    CALL cpu_log( log_point_s(32), 'mpi_alltoall', 'stop' )
1260#else
1261!
1262!-- Next line required to avoid compile error about unused dummy argument in serial mode
1263    i = SIZE( f_in )
1264#endif
1265
1266!
1267!-- Carry out the FFT along x, where all data are present due to the 1d-decomposition along y.
1268!-- Resort the data in a way that y becomes the first index.
1269    CALL cpu_log( log_point_s(4), 'fft_x_1d', 'continue' )
1270
1271    IF ( loop_optimization == 'vector' )  THEN
1272!
1273!--    Code optimized for vector processors
1274!$OMP  PARALLEL PRIVATE ( i, j, k )
1275!$OMP  DO
1276       DO  j = nys, nyn
1277
1278          DO  k = 1, nz
1279             DO  i = 0, nx
1280                work_fftx(i,k,j) = work(j,k,i)
1281             ENDDO
1282          ENDDO
1283
1284          CALL fft_x_m( work_fftx(:,:,j), 'backward' )
1285
1286       ENDDO
1287
1288!$OMP  DO
1289       DO  i = 0, nx
1290          DO  j = nys, nyn
1291             DO  k = 1, nz
1292                f_out(k,j,i) = work_fftx(i,k,j)
1293             ENDDO
1294          ENDDO
1295       ENDDO
1296!$OMP  END PARALLEL
1297
1298    ELSE
1299
1300!
1301!--    Cache optimized code (there might be still a potential for better optimization).
1302!$OMP  PARALLEL PRIVATE (i,j,k)
1303!$OMP  DO
1304       DO  j = nys, nyn
1305          DO  k = 1, nz
1306
1307             DO  i = 0, nx
1308                work_fftx(i,k,j) = work(j,k,i)
1309             ENDDO
1310
1311             CALL fft_x_1d( work_fftx(0:nx,k,j), 'backward' )
1312
1313          ENDDO
1314       ENDDO
1315
1316!$OMP  DO
1317       DO  i = 0, nx
1318          DO  j = nys, nyn
1319             DO  k = 1, nz
1320                f_out(k,j,i) = work_fftx(i,k,j)
1321             ENDDO
1322          ENDDO
1323       ENDDO
1324!$OMP  END PARALLEL
1325
1326    ENDIF
1327    CALL cpu_log( log_point_s(4), 'fft_x_1d', 'stop' )
1328
1329 END SUBROUTINE tr_yx_fftx
1330
1331
1332!--------------------------------------------------------------------------------------------------!
1333! Description:
1334! ------------
1335!> FFT along y, solution of the tridiagonal system and backward FFT for a 1d-decomposition along y.
1336!>
1337!> @warning This subroutine may still not work for hybrid parallelization with OpenMP (for possible
1338!>          necessary changes see the original routine poisfft_hybrid, developed by Klaus Ketelsen,
1339!>          May 2002)
1340!--------------------------------------------------------------------------------------------------!
1341 SUBROUTINE ffty_tri_ffty( ar )
1342
1343
1344    USE control_parameters,                                                                        &
1345        ONLY:  loop_optimization
1346
1347    USE cpulog,                                                                                    &
1348        ONLY:  cpu_log,                                                                            &
1349               log_point_s
1350
1351    USE grid_variables,                                                                            &
1352        ONLY:  ddx2,                                                                               &
1353               ddy2
1354
1355    USE kinds
1356
1357    USE pegrid
1358
1359    IMPLICIT NONE
1360
1361    INTEGER(iwp) ::  i                   !<
1362    INTEGER(iwp) ::  j                   !<
1363    INTEGER(iwp) ::  k                   !<
1364    INTEGER(iwp) ::  m                   !<
1365    INTEGER(iwp) ::  n                   !<
1366!$  INTEGER(iwp) ::  omp_get_thread_num  !<
1367    INTEGER(iwp) ::  tn                  !<
1368
1369    REAL(wp), DIMENSION(0:ny)                          ::  work_ffty  !<
1370    REAL(wp), DIMENSION(0:ny,1:nz)                     ::  work_triy  !<
1371    REAL(wp), DIMENSION(nny,1:nz,nxl_y:nxr_y,pdims(2)) ::  ar         !<
1372    REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE          ::  tri        !<
1373
1374
1375    CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'start' )
1376
1377    ALLOCATE( tri(5,0:ny,0:nz-1,0:threads_per_task-1) )
1378
1379    tn = 0           ! Default thread number in case of one thread
1380!$OMP PARALLEL DO PRIVATE ( i, j, k, m, n, tn, work_ffty, work_triy )
1381    DO  i = nxl_y, nxr_y
1382
1383!$     tn = omp_get_thread_num()
1384
1385       IF ( loop_optimization == 'vector' )  THEN
1386!
1387!--       Code optimized for vector processors
1388          DO  k = 1, nz
1389
1390             m = 0
1391             DO  n = 1, pdims(2)
1392                DO  j = 1, nny
1393                   work_triy(m,k) = ar(j,k,i,n)
1394                   m = m + 1
1395                ENDDO
1396             ENDDO
1397
1398          ENDDO
1399
1400          CALL fft_y_m( work_triy, ny, 'forward' )
1401
1402       ELSE
1403!
1404!--       Cache optimized code
1405          DO  k = 1, nz
1406
1407             m = 0
1408             DO  n = 1, pdims(2)
1409                DO  j = 1, nny
1410                   work_ffty(m) = ar(j,k,i,n)
1411                   m = m + 1
1412                ENDDO
1413             ENDDO
1414
1415             CALL fft_y_1d( work_ffty, 'forward' )
1416
1417             DO  j = 0, ny
1418                work_triy(j,k) = work_ffty(j)
1419             ENDDO
1420
1421          ENDDO
1422
1423       ENDIF
1424
1425!
1426!--    Solve the linear equation system
1427       CALL tridia_1dd( ddy2, ddx2, ny, nx, i, work_triy, tri(:,:,:,tn) )
1428
1429       IF ( loop_optimization == 'vector' )  THEN
1430!
1431!--       Code optimized for vector processors
1432          CALL fft_y_m( work_triy, ny, 'backward' )
1433
1434          DO  k = 1, nz
1435
1436             m = 0
1437             DO  n = 1, pdims(2)
1438                DO  j = 1, nny
1439                   ar(j,k,i,n) = work_triy(m,k)
1440                   m = m + 1
1441                ENDDO
1442             ENDDO
1443
1444          ENDDO
1445
1446       ELSE
1447!
1448!--       Cache optimized code
1449          DO  k = 1, nz
1450
1451             DO  j = 0, ny
1452                work_ffty(j) = work_triy(j,k)
1453             ENDDO
1454
1455             CALL fft_y_1d( work_ffty, 'backward' )
1456
1457             m = 0
1458             DO  n = 1, pdims(2)
1459                DO  j = 1, nny
1460                   ar(j,k,i,n) = work_ffty(m)
1461                   m = m + 1
1462                ENDDO
1463             ENDDO
1464
1465          ENDDO
1466
1467       ENDIF
1468
1469    ENDDO
1470
1471    DEALLOCATE( tri )
1472
1473    CALL cpu_log( log_point_s(39), 'fft_y_1d + tridia', 'stop' )
1474
1475 END SUBROUTINE ffty_tri_ffty
1476
1477 END MODULE poisfft_mod
Note: See TracBrowser for help on using the repository browser.