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

Last change on this file since 684 was 684, checked in by raasch, 10 years ago

last commit documented

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