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

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