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

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

Last commit documented

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