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

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

leapfrog timestep scheme and upstream-spline advection scheme completely removed from the code,
reading of dt_fixed from restart file removed

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