source: palm/trunk/SOURCE/init_pegrid.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 10 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:


Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:


localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

  • Property svn:keywords set to Id
File size: 40.8 KB
Line 
1 SUBROUTINE init_pegrid
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! bc_lr/ns replaced by bc_lr/ns_cyc/dirrad/raddir
7!
8! ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
9! TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values
10!
11! Former revisions:
12! -----------------
13! $Id: init_pegrid.f90 707 2011-03-29 11:39:40Z raasch $
14!
15! 667 2010-12-23 12:06:00Z suehring/gryschka
16! Moved determination of target_id's from init_coupling
17! Determination of parameters needed for coupling (coupling_topology, ngp_a,
18! ngp_o) with different grid/processor-topology in ocean and atmosphere
19! Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
20! The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
21! maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
22! This distinction is due to reasons of data exchange and performance for the
23! normal grid and grids in poismg.
24! The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
25! New MPI-Vectors for data exchange between left and right boundaries added.
26! This is due to reasons of performance (10% faster).
27!
28! 646 2010-12-15 13:03:52Z raasch
29! lctit is now using a 2d decomposition by default
30!
31! 622 2010-12-10 08:08:13Z raasch
32! optional barriers included in order to speed up collective operations
33!
34! 438 2010-02-01 04:32:43Z raasch
35! 2d-decomposition is default for Cray-XT machines
36!
37! 274 2009-03-26 15:11:21Z heinze
38! Output of messages replaced by message handling routine.
39!
40! 206 2008-10-13 14:59:11Z raasch
41! Implementation of a MPI-1 coupling: added __parallel within the __mpi2 part
42! 2d-decomposition is default on SGI-ICE systems
43!
44! 197 2008-09-16 15:29:03Z raasch
45! multigrid levels are limited by subdomains if mg_switch_to_pe0_level = -1,
46! nz is used instead nnz for calculating mg-levels
47! Collect on PE0 horizontal index bounds from all other PEs,
48! broadcast the id of the inflow PE (using the respective communicator)
49!
50! 114 2007-10-10 00:03:15Z raasch
51! Allocation of wall flag arrays for multigrid solver
52!
53! 108 2007-08-24 15:10:38Z letzel
54! Intercommunicator (comm_inter) and derived data type (type_xy) for
55! coupled model runs created, assign coupling_mode_remote,
56! indices nxlu and nysv are calculated (needed for non-cyclic boundary
57! conditions)
58!
59! 82 2007-04-16 15:40:52Z raasch
60! Cpp-directive lcmuk changed to intel_openmp_bug, setting of host on lcmuk by
61! cpp-directive removed
62!
63! 75 2007-03-22 09:54:05Z raasch
64! uxrp, vynp eliminated,
65! dirichlet/neumann changed to dirichlet/radiation, etc.,
66! poisfft_init is only called if fft-solver is switched on
67!
68! RCS Log replace by Id keyword, revision history cleaned up
69!
70! Revision 1.28  2006/04/26 13:23:32  raasch
71! lcmuk does not understand the !$ comment so a cpp-directive is required
72!
73! Revision 1.1  1997/07/24 11:15:09  raasch
74! Initial revision
75!
76!
77! Description:
78! ------------
79! Determination of the virtual processor topology (if not prescribed by the
80! user)and computation of the grid point number and array bounds of the local
81! domains.
82!------------------------------------------------------------------------------!
83
84    USE control_parameters
85    USE fft_xy
86    USE grid_variables
87    USE indices
88    USE pegrid
89    USE poisfft_mod
90    USE poisfft_hybrid_mod
91    USE statistics
92    USE transpose_indices
93
94
95
96    IMPLICIT NONE
97
98    INTEGER ::  gathered_size, i, id_inflow_l, id_recycling_l, ind(5), j, k, &
99                maximum_grid_level_l, mg_switch_to_pe0_level_l, mg_levels_x, &
100                mg_levels_y, mg_levels_z, nnx_y, nnx_z, nny_x, nny_z, nnz_x, &
101                nnz_y, numproc_sqr, nx_total, nxl_l, nxr_l, nyn_l, nys_l,    &
102                nzb_l, nzt_l, omp_get_num_threads, subdomain_size
103
104    INTEGER, DIMENSION(:), ALLOCATABLE ::  ind_all, nxlf, nxrf, nynf, nysf
105
106    INTEGER, DIMENSION(2) :: pdims_remote
107
108    LOGICAL ::  found
109
110!
111!-- Get the number of OpenMP threads
112    !$OMP PARALLEL
113#if defined( __intel_openmp_bug )
114    threads_per_task = omp_get_num_threads()
115#else
116!$  threads_per_task = omp_get_num_threads()
117#endif
118    !$OMP END PARALLEL
119
120
121#if defined( __parallel )
122
123!
124!-- Determine the processor topology or check it, if prescribed by the user
125    IF ( npex == -1  .AND.  npey == -1 )  THEN
126
127!
128!--    Automatic determination of the topology
129!--    The default on SMP- and cluster-hosts is a 1d-decomposition along x
130       IF ( host(1:3) == 'ibm'  .OR.  host(1:3) == 'nec'      .OR. &
131            ( host(1:2) == 'lc'  .AND.  host(3:5) /= 'sgi'  .AND.  &
132              host(3:4) /= 'xt'  .AND.  host(3:5) /= 'tit' )  .OR. &
133             host(1:3) == 'dec' )  THEN
134
135          pdims(1) = numprocs
136          pdims(2) = 1
137
138       ELSE
139
140          numproc_sqr = SQRT( REAL( numprocs ) )
141          pdims(1)    = MAX( numproc_sqr , 1 )
142          DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
143             pdims(1) = pdims(1) - 1
144          ENDDO
145          pdims(2) = numprocs / pdims(1)
146
147       ENDIF
148
149    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
150
151!
152!--    Prescribed by user. Number of processors on the prescribed topology
153!--    must be equal to the number of PEs available to the job
154       IF ( ( npex * npey ) /= numprocs )  THEN
155          WRITE( message_string, * ) 'number of PEs of the prescribed ',      & 
156                 'topology (', npex*npey,') does not match & the number of ', & 
157                 'PEs available to the job (', numprocs, ')'
158          CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
159       ENDIF
160       pdims(1) = npex
161       pdims(2) = npey
162
163    ELSE
164!
165!--    If the processor topology is prescribed by the user, the number of
166!--    PEs must be given in both directions
167       message_string = 'if the processor topology is prescribed by the, ' //  &
168                   ' user& both values of "npex" and "npey" must be given ' // &
169                   'in the &NAMELIST-parameter file'
170       CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
171
172    ENDIF
173
174!
175!-- The hybrid solver can only be used in case of a 1d-decomposition along x
176    IF ( pdims(2) /= 1  .AND.  psolver == 'poisfft_hybrid' )  THEN
177       message_string = 'psolver = "poisfft_hybrid" can only be' // &
178                        '& used in case of a 1d-decomposition along x'
179       CALL message( 'init_pegrid', 'PA0223', 1, 2, 0, 6, 0 )
180    ENDIF
181
182!
183!-- For communication speedup, set barriers in front of collective
184!-- communications by default on SGI-type systems
185    IF ( host(3:5) == 'sgi' )  collective_wait = .TRUE.
186
187!
188!-- If necessary, set horizontal boundary conditions to non-cyclic
189    IF ( .NOT. bc_lr_cyc )  cyclic(1) = .FALSE.
190    IF ( .NOT. bc_ns_cyc )  cyclic(2) = .FALSE.
191
192!
193!-- Create the virtual processor grid
194    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
195                          comm2d, ierr )
196    CALL MPI_COMM_RANK( comm2d, myid, ierr )
197    WRITE (myid_char,'(''_'',I4.4)')  myid
198
199    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
200    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
201    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
202
203!
204!-- Determine sub-topologies for transpositions
205!-- Transposition from z to x:
206    remain_dims(1) = .TRUE.
207    remain_dims(2) = .FALSE.
208    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
209    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
210!
211!-- Transposition from x to y
212    remain_dims(1) = .FALSE.
213    remain_dims(2) = .TRUE.
214    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
215    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
216
217
218!
219!-- Find a grid (used for array d) which will match the transposition demands
220    IF ( grid_matching == 'strict' )  THEN
221
222       nxa = nx;  nya = ny;  nza = nz
223
224    ELSE
225
226       found = .FALSE.
227   xn: DO  nxa = nx, 2*nx
228!
229!--       Meet conditions for nx
230          IF ( MOD( nxa+1, pdims(1) ) /= 0 .OR. &
231               MOD( nxa+1, pdims(2) ) /= 0 )  CYCLE xn
232
233      yn: DO  nya = ny, 2*ny
234!
235!--          Meet conditions for ny
236             IF ( MOD( nya+1, pdims(2) ) /= 0 .OR. &
237                  MOD( nya+1, pdims(1) ) /= 0 )  CYCLE yn
238
239
240         zn: DO  nza = nz, 2*nz
241!
242!--             Meet conditions for nz
243                IF ( ( MOD( nza, pdims(1) ) /= 0  .AND.  pdims(1) /= 1  .AND. &
244                       pdims(2) /= 1 )  .OR.                                  &
245                     ( MOD( nza, pdims(2) ) /= 0  .AND.  dt_dosp /= 9999999.9 &
246                     ) )  THEN
247                   CYCLE zn
248                ELSE
249                   found = .TRUE.
250                   EXIT xn
251                ENDIF
252
253             ENDDO zn
254
255          ENDDO yn
256
257       ENDDO xn
258
259       IF ( .NOT. found )  THEN
260          message_string = 'no matching grid for transpositions found'
261          CALL message( 'init_pegrid', 'PA0224', 1, 2, 0, 6, 0 )
262       ENDIF
263
264    ENDIF
265
266!
267!-- Calculate array bounds in x-direction for every PE.
268!-- The last PE along x may get less grid points than the others
269    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1), &
270              nysf(0:pdims(2)-1), nnx_pe(0:pdims(1)-1), nny_pe(0:pdims(2)-1) )
271
272    IF ( MOD( nxa+1 , pdims(1) ) /= 0 )  THEN
273       WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ',&
274                               'is not an& integral divisor of the number ',  &
275                               'processors (', pdims(1),')'
276       CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
277    ELSE
278       nnx  = ( nxa + 1 ) / pdims(1)
279       IF ( nnx*pdims(1) - ( nx + 1) > nnx )  THEN
280          WRITE( message_string, * ) 'x-direction: nx does not match the',    & 
281                       'requirements given by the number of PEs &used',       &
282                       '& please use nx = ', nx - ( pdims(1) - ( nnx*pdims(1) &
283                                   - ( nx + 1 ) ) ), ' instead of nx =', nx
284          CALL message( 'init_pegrid', 'PA0226', 1, 2, 0, 6, 0 )
285       ENDIF
286    ENDIF   
287
288!
289!-- Left and right array bounds, number of gridpoints
290    DO  i = 0, pdims(1)-1
291       nxlf(i)   = i * nnx
292       nxrf(i)   = ( i + 1 ) * nnx - 1
293       nnx_pe(i) = MIN( nx, nxrf(i) ) - nxlf(i) + 1
294    ENDDO
295
296!
297!-- Calculate array bounds in y-direction for every PE.
298    IF ( MOD( nya+1 , pdims(2) ) /= 0 )  THEN
299       WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
300                           'is not an& integral divisor of the number of',     &
301                           'processors (', pdims(2),')'
302       CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
303    ELSE
304       nny  = ( nya + 1 ) / pdims(2)
305       IF ( nny*pdims(2) - ( ny + 1) > nny )  THEN
306          WRITE( message_string, * ) 'y-direction: ny does not match the',    &
307                       'requirements given by the number of PEs &used ',      &
308                       '& please use ny = ', ny - ( pdims(2) - ( nnx*pdims(2) &
309                                     - ( ny + 1 ) ) ), ' instead of ny =', ny
310          CALL message( 'init_pegrid', 'PA0228', 1, 2, 0, 6, 0 )
311       ENDIF
312    ENDIF   
313
314!
315!-- South and north array bounds
316    DO  j = 0, pdims(2)-1
317       nysf(j)   = j * nny
318       nynf(j)   = ( j + 1 ) * nny - 1
319       nny_pe(j) = MIN( ny, nynf(j) ) - nysf(j) + 1
320    ENDDO
321
322!
323!-- Local array bounds of the respective PEs
324    nxl  = nxlf(pcoord(1))
325    nxra = nxrf(pcoord(1))
326    nxr  = MIN( nx, nxra )
327    nys  = nysf(pcoord(2))
328    nyna = nynf(pcoord(2))
329    nyn  = MIN( ny, nyna )
330    nzb  = 0
331    nzta = nza
332    nzt  = MIN( nz, nzta )
333    nnz  = nza
334
335!
336!-- Set switches to define if the PE is situated at the border of the virtual
337!-- processor grid
338    IF ( nxl == 0 )   left_border_pe  = .TRUE.
339    IF ( nxr == nx )  right_border_pe = .TRUE.
340    IF ( nys == 0 )   south_border_pe = .TRUE.
341    IF ( nyn == ny )  north_border_pe = .TRUE.
342
343!
344!-- Calculate array bounds and gridpoint numbers for the transposed arrays
345!-- (needed in the pressure solver)
346!-- For the transposed arrays, cyclic boundaries as well as top and bottom
347!-- boundaries are omitted, because they are obstructive to the transposition
348
349!
350!-- 1. transposition  z --> x
351!-- This transposition is not neccessary in case of a 1d-decomposition along x,
352!-- except that the uptream-spline method is switched on
353    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
354         scalar_advec == 'ups-scheme' )  THEN
355
356       IF ( pdims(2) == 1  .AND. ( momentum_advec == 'ups-scheme'  .OR. &
357            scalar_advec == 'ups-scheme' ) )  THEN
358          message_string = '1d-decomposition along x ' // &
359                           'chosen but nz restrictions may occur' // &
360                           '& since ups-scheme is activated'
361          CALL message( 'init_pegrid', 'PA0229', 0, 1, 0, 6, 0 )
362       ENDIF
363       nys_x  = nys
364       nyn_xa = nyna
365       nyn_x  = nyn
366       nny_x  = nny
367       IF ( MOD( nza , pdims(1) ) /= 0 )  THEN
368          WRITE( message_string, * ) 'transposition z --> x:',                &
369                       '&nz=',nz,' is not an integral divisior of pdims(1)=', &
370                                                                   pdims(1)
371          CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
372       ENDIF
373       nnz_x  = nza / pdims(1)
374       nzb_x  = 1 + myidx * nnz_x
375       nzt_xa = ( myidx + 1 ) * nnz_x
376       nzt_x  = MIN( nzt, nzt_xa )
377
378       sendrecvcount_zx = nnx * nny * nnz_x
379
380    ELSE
381!
382!---   Setting of dummy values because otherwise variables are undefined in
383!---   the next step  x --> y
384!---   WARNING: This case has still to be clarified!!!!!!!!!!!!
385       nnz_x  = 1
386       nzb_x  = 1
387       nzt_xa = 1
388       nzt_x  = 1
389       nny_x  = nny
390
391    ENDIF
392
393!
394!-- 2. transposition  x --> y
395    nnz_y  = nnz_x
396    nzb_y  = nzb_x
397    nzt_ya = nzt_xa
398    nzt_y  = nzt_x
399    IF ( MOD( nxa+1 , pdims(2) ) /= 0 )  THEN
400       WRITE( message_string, * ) 'transposition x --> y:',                &
401                         '&nx+1=',nx+1,' is not an integral divisor of ',&
402                         'pdims(2)=',pdims(2)
403       CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
404    ENDIF
405    nnx_y = (nxa+1) / pdims(2)
406    nxl_y = myidy * nnx_y
407    nxr_ya = ( myidy + 1 ) * nnx_y - 1
408    nxr_y  = MIN( nx, nxr_ya )
409
410    sendrecvcount_xy = nnx_y * nny_x * nnz_y
411
412!
413!-- 3. transposition  y --> z  (ELSE:  x --> y  in case of 1D-decomposition
414!-- along x)
415    IF ( pdims(2) /= 1  .OR.  momentum_advec == 'ups-scheme'  .OR. &
416         scalar_advec == 'ups-scheme' )  THEN
417!
418!--    y --> z
419!--    This transposition is not neccessary in case of a 1d-decomposition
420!--    along x, except that the uptream-spline method is switched on
421       nnx_z  = nnx_y
422       nxl_z  = nxl_y
423       nxr_za = nxr_ya
424       nxr_z  = nxr_y
425       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
426          WRITE( message_string, * ) 'transposition y --> z:',            &
427                            '& ny+1=',ny+1,' is not an integral divisor of ',&
428                            'pdims(1)=',pdims(1)
429          CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
430       ENDIF
431       nny_z  = (nya+1) / pdims(1)
432       nys_z  = myidx * nny_z
433       nyn_za = ( myidx + 1 ) * nny_z - 1
434       nyn_z  = MIN( ny, nyn_za )
435
436       sendrecvcount_yz = nnx_y * nny_z * nnz_y
437
438    ELSE
439!
440!--    x --> y. This condition must be fulfilled for a 1D-decomposition along x
441       IF ( MOD( nya+1 , pdims(1) ) /= 0 )  THEN
442          WRITE( message_string, * ) 'transposition x --> y:',               &
443                            '& ny+1=',ny+1,' is not an integral divisor of ',&
444                            'pdims(1)=',pdims(1)
445          CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
446       ENDIF
447
448    ENDIF
449
450!
451!-- Indices for direct transpositions z --> y (used for calculating spectra)
452    IF ( dt_dosp /= 9999999.9 )  THEN
453       IF ( MOD( nza, pdims(2) ) /= 0 )  THEN
454          WRITE( message_string, * ) 'direct transposition z --> y (needed ', &
455                    'for spectra):& nz=',nz,' is not an integral divisor of ',&
456                    'pdims(2)=',pdims(2)
457          CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
458       ELSE
459          nxl_yd  = nxl
460          nxr_yda = nxra
461          nxr_yd  = nxr
462          nzb_yd  = 1 + myidy * ( nza / pdims(2) )
463          nzt_yda = ( myidy + 1 ) * ( nza / pdims(2) )
464          nzt_yd  = MIN( nzt, nzt_yda )
465
466          sendrecvcount_zyd = nnx * nny * ( nza / pdims(2) )
467       ENDIF
468    ENDIF
469
470!
471!-- Indices for direct transpositions y --> x (they are only possible in case
472!-- of a 1d-decomposition along x)
473    IF ( pdims(2) == 1 )  THEN
474       nny_x  = nny / pdims(1)
475       nys_x  = myid * nny_x
476       nyn_xa = ( myid + 1 ) * nny_x - 1
477       nyn_x  = MIN( ny, nyn_xa )
478       nzb_x  = 1
479       nzt_xa = nza
480       nzt_x  = nz
481       sendrecvcount_xy = nnx * nny_x * nza
482    ENDIF
483
484!
485!-- Indices for direct transpositions x --> y (they are only possible in case
486!-- of a 1d-decomposition along y)
487    IF ( pdims(1) == 1 )  THEN
488       nnx_y  = nnx / pdims(2)
489       nxl_y  = myid * nnx_y
490       nxr_ya = ( myid + 1 ) * nnx_y - 1
491       nxr_y  = MIN( nx, nxr_ya )
492       nzb_y  = 1
493       nzt_ya = nza
494       nzt_y  = nz
495       sendrecvcount_xy = nnx_y * nny * nza
496    ENDIF
497
498!
499!-- Arrays for storing the array bounds are needed any more
500    DEALLOCATE( nxlf , nxrf , nynf , nysf )
501
502!
503!-- Collect index bounds from other PEs (to be written to restart file later)
504    ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
505
506    IF ( myid == 0 )  THEN
507
508       hor_index_bounds(1,0) = nxl
509       hor_index_bounds(2,0) = nxr
510       hor_index_bounds(3,0) = nys
511       hor_index_bounds(4,0) = nyn
512
513!
514!--    Receive data from all other PEs
515       DO  i = 1, numprocs-1
516          CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
517                         ierr )
518          hor_index_bounds(:,i) = ibuf(1:4)
519       ENDDO
520
521    ELSE
522!
523!--    Send index bounds to PE0
524       ibuf(1) = nxl
525       ibuf(2) = nxr
526       ibuf(3) = nys
527       ibuf(4) = nyn
528       CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
529
530    ENDIF
531
532#if defined( __print )
533!
534!-- Control output
535    IF ( myid == 0 )  THEN
536       PRINT*, '*** processor topology ***'
537       PRINT*, ' '
538       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
539               &'   nys: nyn'
540       PRINT*, '------------------------------------------------------------',&
541               &'-----------'
542       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
543                       myidx, myidy, nxl, nxr, nys, nyn
5441000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
545               2(2X,I4,':',I4))
546
547!
548!--    Receive data from the other PEs
549       DO  i = 1,numprocs-1
550          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
551                         ierr )
552          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
553       ENDDO
554    ELSE
555
556!
557!--    Send data to PE0
558       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
559       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
560       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
561       ibuf(12) = nyn
562       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
563    ENDIF
564#endif
565
566#if defined( __parallel )
567#if defined( __mpi2 )
568!
569!-- In case of coupled runs, get the port name on PE0 of the atmosphere model
570!-- and pass it to PE0 of the ocean model
571    IF ( myid == 0 )  THEN
572
573       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
574
575          CALL MPI_OPEN_PORT( MPI_INFO_NULL, port_name, ierr )
576!
577!--       TEST OUTPUT (TO BE REMOVED)
578          WRITE(9,*)  TRIM( coupling_mode ),  &
579               ', ierr after MPI_OPEN_PORT: ', ierr
580          CALL LOCAL_FLUSH( 9 )
581
582          CALL MPI_PUBLISH_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, &
583                                 ierr )
584!
585!--       TEST OUTPUT (TO BE REMOVED)
586          WRITE(9,*)  TRIM( coupling_mode ),  &
587               ', ierr after MPI_PUBLISH_NAME: ', ierr
588          CALL LOCAL_FLUSH( 9 )
589
590!
591!--       Write a flag file for the ocean model and the other atmosphere
592!--       processes.
593!--       There seems to be a bug in MPICH2 which causes hanging processes
594!--       in case that execution of LOOKUP_NAME is continued too early
595!--       (i.e. before the port has been created)
596          OPEN( 90, FILE='COUPLING_PORT_OPENED', FORM='FORMATTED' )
597          WRITE ( 90, '(''TRUE'')' )
598          CLOSE ( 90 )
599
600       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
601
602!
603!--       Continue only if the atmosphere model has created the port.
604!--       There seems to be a bug in MPICH2 which causes hanging processes
605!--       in case that execution of LOOKUP_NAME is continued too early
606!--       (i.e. before the port has been created)
607          INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
608          DO WHILE ( .NOT. found )
609             INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
610          ENDDO
611
612          CALL MPI_LOOKUP_NAME( 'palm_coupler', MPI_INFO_NULL, port_name, ierr )
613!
614!--       TEST OUTPUT (TO BE REMOVED)
615          WRITE(9,*)  TRIM( coupling_mode ),  &
616               ', ierr after MPI_LOOKUP_NAME: ', ierr
617          CALL LOCAL_FLUSH( 9 )
618
619
620       ENDIF
621
622    ENDIF
623
624!
625!-- In case of coupled runs, establish the connection between the atmosphere
626!-- and the ocean model and define the intercommunicator (comm_inter)
627    CALL MPI_BARRIER( comm2d, ierr )
628    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
629
630       PRINT*, '... before COMM_ACCEPT'
631       CALL MPI_COMM_ACCEPT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
632                             comm_inter, ierr )
633       PRINT*, '--- ierr = ', ierr
634       PRINT*, '--- comm_inter atmosphere = ', comm_inter
635
636       coupling_mode_remote = 'ocean_to_atmosphere'
637
638    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
639
640       IF ( myid == 0 )  PRINT*, '*** read: ', port_name, '  ierr = ', ierr
641       PRINT*, '... before COMM_CONNECT'
642       CALL MPI_COMM_CONNECT( port_name, MPI_INFO_NULL, 0, MPI_COMM_WORLD, &
643                              comm_inter, ierr )
644       PRINT*, '--- ierr = ', ierr
645       PRINT*, '--- comm_inter ocean      = ', comm_inter
646
647       coupling_mode_remote = 'atmosphere_to_ocean'
648
649    ENDIF
650#endif
651
652!
653!-- Determine the number of ghost points
654    IF (scalar_advec == 'ws-scheme' .OR. momentum_advec == 'ws-scheme') THEN
655       nbgp = 3
656    ELSE
657       nbgp = 1
658    END IF 
659
660!
661!-- In case of coupled runs, create a new MPI derived datatype for the
662!-- exchange of surface (xy) data .
663!-- Gridpoint number for the exchange of ghost points (xy-plane)
664
665    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
666
667!
668!-- Define a new MPI derived datatype for the exchange of ghost points in
669!-- y-direction for 2D-arrays (line)
670    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
671    CALL MPI_TYPE_COMMIT( type_xy, ierr )
672
673
674    IF ( TRIM( coupling_mode ) .NE. 'uncoupled' ) THEN
675   
676!
677!--    Pass the number of grid points of the atmosphere model to
678!--    the ocean model and vice versa
679       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
680
681          nx_a = nx
682          ny_a = ny
683
684          IF ( myid == 0 ) THEN
685             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, &
686                            comm_inter, ierr )
687             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, &
688                            comm_inter, ierr )
689             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, &
690                            comm_inter, ierr )
691             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, &
692                            comm_inter, status, ierr )
693             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, &
694                            comm_inter, status, ierr )
695             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6, &
696                            comm_inter, status, ierr )
697          ENDIF
698
699          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr)
700          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr) 
701          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr)
702       
703       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
704
705          nx_o = nx
706          ny_o = ny
707
708          IF ( myid == 0 ) THEN
709             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, &
710                            comm_inter, status, ierr )
711             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, &
712                            comm_inter, status, ierr )
713             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, &
714                            comm_inter, status, ierr )
715             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, &
716                            comm_inter, ierr )
717             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, &
718                            comm_inter, ierr )
719             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, &
720                            comm_inter, ierr )
721          ENDIF
722
723          CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
724          CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr) 
725          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr) 
726
727       ENDIF
728 
729       ngp_a = (nx_a+1+2*nbgp)*(ny_a+1+2*nbgp)
730       ngp_o = (nx_o+1+2*nbgp)*(ny_o+1+2*nbgp)
731
732!
733!--    determine if the horizontal grid and the number of PEs
734!--    in ocean and atmosphere is same or not
735!--    (different number of PEs still not implemented)
736       IF ( nx_o == nx_a .AND. ny_o == ny_a .AND.  &
737            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
738       THEN
739          coupling_topology = 0
740       ELSE
741          coupling_topology = 1
742       ENDIF
743
744!
745!--    Determine the target PEs for the exchange between ocean and
746!--    atmosphere (comm2d)
747       IF ( coupling_topology == 0) THEN
748          IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
749             target_id = myid + numprocs
750          ELSE
751             target_id = myid
752          ENDIF
753
754       ELSE
755
756!
757!--       In case of nonequivalent topology in ocean and atmosphere only for
758!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
759!--       data echxchange between ocean and atmosphere will be done only by
760!--       those PEs.   
761          IF ( myid == 0 ) THEN
762             IF ( TRIM( coupling_mode ) .EQ. 'atmosphere_to_ocean' ) THEN
763                target_id = numprocs
764             ELSE
765                target_id = 0
766             ENDIF
767 print*, coupling_mode, myid, " -> ", target_id, "numprocs: ", numprocs
768          ENDIF
769       ENDIF
770
771    ENDIF
772
773
774#endif
775
776#else
777
778!
779!-- Array bounds when running on a single PE (respectively a non-parallel
780!-- machine)
781    nxl  = 0
782    nxr  = nx
783    nxra = nx
784    nnx  = nxr - nxl + 1
785    nys  = 0
786    nyn  = ny
787    nyna = ny
788    nny  = nyn - nys + 1
789    nzb  = 0
790    nzt  = nz
791    nzta = nz
792    nnz  = nz
793
794    ALLOCATE( hor_index_bounds(4,0:0) )
795    hor_index_bounds(1,0) = nxl
796    hor_index_bounds(2,0) = nxr
797    hor_index_bounds(3,0) = nys
798    hor_index_bounds(4,0) = nyn
799
800!
801!-- Array bounds for the pressure solver (in the parallel code, these bounds
802!-- are the ones for the transposed arrays)
803    nys_x  = nys
804    nyn_x  = nyn
805    nyn_xa = nyn
806    nzb_x  = nzb + 1
807    nzt_x  = nzt
808    nzt_xa = nzt
809
810    nxl_y  = nxl
811    nxr_y  = nxr
812    nxr_ya = nxr
813    nzb_y  = nzb + 1
814    nzt_y  = nzt
815    nzt_ya = nzt
816
817    nxl_z  = nxl
818    nxr_z  = nxr
819    nxr_za = nxr
820    nys_z  = nys
821    nyn_z  = nyn
822    nyn_za = nyn
823
824#endif
825
826!
827!-- Calculate number of grid levels necessary for the multigrid poisson solver
828!-- as well as the gridpoint indices on each level
829    IF ( psolver == 'multigrid' )  THEN
830
831!
832!--    First calculate number of possible grid levels for the subdomains
833       mg_levels_x = 1
834       mg_levels_y = 1
835       mg_levels_z = 1
836
837       i = nnx
838       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
839          i = i / 2
840          mg_levels_x = mg_levels_x + 1
841       ENDDO
842
843       j = nny
844       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
845          j = j / 2
846          mg_levels_y = mg_levels_y + 1
847       ENDDO
848
849       k = nz    ! do not use nnz because it might be > nz due to transposition
850                 ! requirements
851       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
852          k = k / 2
853          mg_levels_z = mg_levels_z + 1
854       ENDDO
855
856       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
857
858!
859!--    Find out, if the total domain allows more levels. These additional
860!--    levels are processed on PE0 only.
861       IF ( numprocs > 1  .AND.  mg_switch_to_pe0_level /= -1 )  THEN
862          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
863             mg_switch_to_pe0_level_l = maximum_grid_level
864
865             mg_levels_x = 1
866             mg_levels_y = 1
867
868             i = nx+1
869             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
870                i = i / 2
871                mg_levels_x = mg_levels_x + 1
872             ENDDO
873
874             j = ny+1
875             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
876                j = j / 2
877                mg_levels_y = mg_levels_y + 1
878             ENDDO
879
880             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
881
882             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
883                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
884                                           mg_switch_to_pe0_level_l + 1
885             ELSE
886                mg_switch_to_pe0_level_l = 0
887             ENDIF
888          ELSE
889             mg_switch_to_pe0_level_l = 0
890             maximum_grid_level_l = maximum_grid_level
891          ENDIF
892
893!
894!--       Use switch level calculated above only if it is not pre-defined
895!--       by user
896          IF ( mg_switch_to_pe0_level == 0 )  THEN
897
898             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
899                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
900                maximum_grid_level     = maximum_grid_level_l
901             ENDIF
902
903          ELSE
904!
905!--          Check pre-defined value and reset to default, if neccessary
906             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.  &
907                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
908                message_string = 'mg_switch_to_pe0_level ' // &
909                                 'out of range and reset to default (=0)'
910                CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
911                mg_switch_to_pe0_level = 0
912             ELSE
913!
914!--             Use the largest number of possible levels anyway and recalculate
915!--             the switch level to this largest number of possible values
916                maximum_grid_level = maximum_grid_level_l
917
918             ENDIF
919          ENDIF
920
921       ENDIF
922
923       ALLOCATE( grid_level_count(maximum_grid_level),                   &
924                 nxl_mg(maximum_grid_level), nxr_mg(maximum_grid_level), &
925                 nyn_mg(maximum_grid_level), nys_mg(maximum_grid_level), &
926                 nzt_mg(maximum_grid_level) )
927
928       grid_level_count = 0
929       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
930
931       DO  i = maximum_grid_level, 1 , -1
932
933          IF ( i == mg_switch_to_pe0_level )  THEN
934#if defined( __parallel )
935!
936!--          Save the grid size of the subdomain at the switch level, because
937!--          it is needed in poismg.
938!--          Array bounds of the local subdomain grids are gathered on PE0
939             ind(1) = nxl_l; ind(2) = nxr_l
940             ind(3) = nys_l; ind(4) = nyn_l
941             ind(5) = nzt_l
942             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
943             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
944                                 MPI_INTEGER, comm2d, ierr )
945             DO  j = 0, numprocs-1
946                DO  k = 1, 5
947                   mg_loc_ind(k,j) = ind_all(k+j*5)
948                ENDDO
949             ENDDO
950             DEALLOCATE( ind_all )
951!
952!--          Calculate the grid size of the total domain gathered on PE0
953             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
954             nxl_l = 0
955             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
956             nys_l = 0
957!
958!--          The size of this gathered array must not be larger than the
959!--          array tend, which is used in the multigrid scheme as a temporary
960!--          array
961             subdomain_size = ( nxr - nxl + 3 )     * ( nyn - nys + 3 )     * &
962                              ( nzt - nzb + 2 )
963             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) * &
964                              ( nzt_l - nzb + 2 )
965
966             IF ( gathered_size > subdomain_size )  THEN
967                message_string = 'not enough memory for storing ' // &
968                                 'gathered multigrid data on PE0'
969                CALL message( 'init_pegrid', 'PA0236', 1, 2, 0, 6, 0 )
970             ENDIF
971#else
972             message_string = 'multigrid gather/scatter impossible ' // &
973                          'in non parallel mode'
974             CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
975#endif
976          ENDIF
977
978          nxl_mg(i) = nxl_l
979          nxr_mg(i) = nxr_l
980          nys_mg(i) = nys_l
981          nyn_mg(i) = nyn_l
982          nzt_mg(i) = nzt_l
983
984          nxl_l = nxl_l / 2 
985          nxr_l = nxr_l / 2
986          nys_l = nys_l / 2 
987          nyn_l = nyn_l / 2 
988          nzt_l = nzt_l / 2 
989       ENDDO
990
991    ELSE
992
993       maximum_grid_level = 0
994
995    ENDIF
996
997    grid_level = maximum_grid_level
998
999#if defined( __parallel )
1000!
1001!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
1002    ngp_y  = nyn - nys + 1 + 2 * nbgp
1003
1004!
1005!-- Define a new MPI derived datatype for the exchange of ghost points in
1006!-- y-direction for 2D-arrays (line)
1007    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x, ierr )
1008    CALL MPI_TYPE_COMMIT( type_x, ierr )
1009    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER, type_x_int, ierr )
1010    CALL MPI_TYPE_COMMIT( type_x_int, ierr )
1011
1012    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
1013    CALL MPI_TYPE_COMMIT( type_y, ierr )
1014    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int, ierr )
1015    CALL MPI_TYPE_COMMIT( type_y_int, ierr )
1016
1017
1018!
1019!-- Calculate gridpoint numbers for the exchange of ghost points along x
1020!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
1021!-- exchange of ghost points in y-direction (xz-plane).
1022!-- Do these calculations for the model grid and (if necessary) also
1023!-- for the coarser grid levels used in the multigrid method
1024    ALLOCATE ( ngp_yz(0:maximum_grid_level), type_xz(0:maximum_grid_level),&
1025               type_yz(0:maximum_grid_level) )
1026
1027    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
1028!
1029!-- Discern between the model grid, which needs nbgp ghost points and
1030!-- grid levels for the multigrid scheme. In the latter case only one
1031!-- ghost point is necessary.
1032!-- First definition of mpi-vectors for exchange of ghost layers on normal
1033!-- grid. The following loop is needed for data exchange in poismg.f90.
1034!
1035!-- Determine number of grid points of yz-layer for exchange
1036    ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1037!
1038!-- Define a new mpi datatype for the exchange of left - right boundaries.
1039!-- Indeed the data are connected in the physical memory and no mpi-vector
1040!-- is necessary, but the data exchange between left and right PE's using
1041!-- mpi-vectors is 10% faster than without.
1042    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz(0), &
1043                             MPI_REAL, type_xz(0), ierr )
1044    CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
1045
1046    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), ierr) 
1047    CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
1048!
1049!-- Definition of mpi-vectors for multigrid
1050    IF ( psolver == 'multigrid' )  THEN
1051!   
1052!--   The definition of mpi-vectors as aforementioned, but only 1 ghost point is used.
1053       DO i = maximum_grid_level, 1 , -1
1054          ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1055
1056          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i), &
1057                             MPI_REAL, type_xz(i), ierr )
1058          CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
1059
1060          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), ierr)
1061          CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
1062
1063          nxl_l = nxl_l / 2
1064          nxr_l = nxr_l / 2
1065          nys_l = nys_l / 2
1066          nyn_l = nyn_l / 2
1067          nzt_l = nzt_l / 2
1068       ENDDO
1069    END IF
1070#endif
1071
1072#if defined( __parallel )
1073!
1074!-- Setting of flags for inflow/outflow conditions in case of non-cyclic
1075!-- horizontal boundary conditions.
1076    IF ( pleft == MPI_PROC_NULL )  THEN
1077       IF ( bc_lr_dirrad )  THEN
1078          inflow_l  = .TRUE.
1079       ELSEIF ( bc_lr_raddir )  THEN
1080          outflow_l = .TRUE.
1081       ENDIF
1082    ENDIF
1083
1084    IF ( pright == MPI_PROC_NULL )  THEN
1085       IF ( bc_lr_dirrad )  THEN
1086          outflow_r = .TRUE.
1087       ELSEIF ( bc_lr_raddir )  THEN
1088          inflow_r  = .TRUE.
1089       ENDIF
1090    ENDIF
1091
1092    IF ( psouth == MPI_PROC_NULL )  THEN
1093       IF ( bc_ns_dirrad )  THEN
1094          outflow_s = .TRUE.
1095       ELSEIF ( bc_ns_raddir )  THEN
1096          inflow_s  = .TRUE.
1097       ENDIF
1098    ENDIF
1099
1100    IF ( pnorth == MPI_PROC_NULL )  THEN
1101       IF ( bc_ns_dirrad )  THEN
1102          inflow_n  = .TRUE.
1103       ELSEIF ( bc_ns_raddir )  THEN
1104          outflow_n = .TRUE.
1105       ENDIF
1106    ENDIF
1107
1108!
1109!-- Broadcast the id of the inflow PE
1110    IF ( inflow_l )  THEN
1111       id_inflow_l = myidx
1112    ELSE
1113       id_inflow_l = 0
1114    ENDIF
1115    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1116    CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
1117                        comm1dx, ierr )
1118
1119!
1120!-- Broadcast the id of the recycling plane
1121!-- WARNING: needs to be adjusted in case of inflows other than from left side!
1122    IF ( ( recycling_width / dx ) >= nxl  .AND. &
1123         ( recycling_width / dx ) <= nxr )  THEN
1124       id_recycling_l = myidx
1125    ELSE
1126       id_recycling_l = 0
1127    ENDIF
1128    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1129    CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
1130                        comm1dx, ierr )
1131
1132#else
1133    IF ( bc_lr_dirrad )  THEN
1134       inflow_l  = .TRUE.
1135       outflow_r = .TRUE.
1136    ELSEIF ( bc_lr_raddir )  THEN
1137       outflow_l = .TRUE.
1138       inflow_r  = .TRUE.
1139    ENDIF
1140
1141    IF ( bc_ns_dirrad )  THEN
1142       inflow_n  = .TRUE.
1143       outflow_s = .TRUE.
1144    ELSEIF ( bc_ns_raddir )  THEN
1145       outflow_n = .TRUE.
1146       inflow_s  = .TRUE.
1147    ENDIF
1148#endif
1149!
1150!-- At the outflow, u or v, respectively, have to be calculated for one more
1151!-- grid point.
1152    IF ( outflow_l )  THEN
1153       nxlu = nxl + 1
1154    ELSE
1155       nxlu = nxl
1156    ENDIF
1157    IF ( outflow_s )  THEN
1158       nysv = nys + 1
1159    ELSE
1160       nysv = nys
1161    ENDIF
1162
1163    IF ( psolver == 'poisfft_hybrid' )  THEN
1164       CALL poisfft_hybrid_ini
1165    ELSEIF ( psolver == 'poisfft' )  THEN
1166       CALL poisfft_init
1167    ENDIF
1168
1169!
1170!-- Allocate wall flag arrays used in the multigrid solver
1171    IF ( psolver == 'multigrid' )  THEN
1172
1173       DO  i = maximum_grid_level, 1, -1
1174
1175           SELECT CASE ( i )
1176
1177              CASE ( 1 )
1178                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
1179                                        nys_mg(i)-1:nyn_mg(i)+1, &
1180                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1181
1182              CASE ( 2 )
1183                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
1184                                        nys_mg(i)-1:nyn_mg(i)+1, &
1185                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1186
1187              CASE ( 3 )
1188                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
1189                                        nys_mg(i)-1:nyn_mg(i)+1, &
1190                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1191
1192              CASE ( 4 )
1193                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
1194                                        nys_mg(i)-1:nyn_mg(i)+1, &
1195                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1196
1197              CASE ( 5 )
1198                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
1199                                        nys_mg(i)-1:nyn_mg(i)+1, &
1200                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1201
1202              CASE ( 6 )
1203                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
1204                                        nys_mg(i)-1:nyn_mg(i)+1, &
1205                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1206
1207              CASE ( 7 )
1208                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
1209                                        nys_mg(i)-1:nyn_mg(i)+1, &
1210                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1211
1212              CASE ( 8 )
1213                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
1214                                        nys_mg(i)-1:nyn_mg(i)+1, &
1215                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1216
1217              CASE ( 9 )
1218                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
1219                                        nys_mg(i)-1:nyn_mg(i)+1, &
1220                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1221
1222              CASE ( 10 )
1223                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
1224                                        nys_mg(i)-1:nyn_mg(i)+1, &
1225                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1226
1227              CASE DEFAULT
1228                 message_string = 'more than 10 multigrid levels'
1229                 CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
1230
1231          END SELECT
1232
1233       ENDDO
1234
1235    ENDIF
1236
1237 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.