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

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

last commit documented

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