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

Last change on this file since 778 was 778, checked in by fricke, 10 years ago

Modifications of the multigrid pressure solver

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