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

Last change on this file since 668 was 668, checked in by suehring, 11 years ago

last commit documented

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