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

Last change on this file since 807 was 807, checked in by maronga, 10 years ago

new utility check_namelist_files implemented

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