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

Last change on this file since 979 was 979, checked in by fricke, 9 years ago

last commit documented

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