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

Last change on this file since 1002 was 1002, checked in by raasch, 9 years ago

last commit documented

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