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

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

Bugfix: cpp directives .NOT., .AND. replaced by !,&&. Minor bugfixes in mrungui

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