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

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

New:
---

Optional barriers included in order to speed up collective operations
MPI_ALLTOALL and MPI_ALLREDUCE. This feature is controlled with new initial
parameter collective_wait. Default is .FALSE, but .TRUE. on SGI-type
systems. (advec_particles, advec_s_bc, buoyancy, check_for_restart,
cpu_statistics, data_output_2d, data_output_ptseries, flow_statistics,
global_min_max, inflow_turbulence, init_3d_model, init_particles, init_pegrid,
init_slope, parin, pres, poismg, set_particle_attributes, timestep,
read_var_list, user_statistics, write_compressed, write_var_list)

Adjustments for Kyushu Univ. (lcrte, ibmku). Concerning hybrid
(MPI/openMP) runs, the number of openMP threads per MPI tasks can now
be given as an argument to mrun-option -O. (mbuild, mrun, subjob)

Changed:


Initialization of the module command changed for SGI-ICE/lcsgi (mbuild, subjob)

Errors:


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