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

Last change on this file since 761 was 761, checked in by suehring, 13 years ago

Bugfix: Avoid divisions by zero in case of using a 'neumann' bc for the pressure at the top of the model domain.

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