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

Last change on this file since 1094 was 1093, checked in by raasch, 11 years ago

last commit documented

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