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

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

last commit documented

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