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

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

New:
---

GPU porting of pres, swap_timelevel. Adjustments of openACC directives.
Further porting of poisfft, which now runs completely on GPU without any
host/device data transfer for serial an parallel runs (but parallel runs
require data transfer before and after the MPI transpositions).
GPU-porting of tridiagonal solver:
tridiagonal routines split into extermal subroutines (instead using CONTAINS),
no distinction between parallel/non-parallel in poisfft and tridia any more,
tridia routines moved to end of file because of probable bug in PGI compiler
(otherwise "invalid device function" is indicated during runtime).
(cuda_fft_interfaces, fft_xy, flow_statistics, init_3d_model, palm, poisfft, pres, prognostic_equations, swap_timelevel, time_integration, transpose)
output of accelerator board information. (header)

optimization of tridia routines: constant elements and coefficients of tri are
stored in seperate arrays ddzuw and tric, last dimension of tri reduced from 5 to 2,
(init_grid, init_3d_model, modules, palm, poisfft)

poisfft_init is now called internally from poisfft,
(Makefile, Makefile_check, init_pegrid, poisfft, poisfft_hybrid)

CPU-time per grid point and timestep is output to CPU_MEASURES file
(cpu_statistics, modules, time_integration)

Changed:


resorting from/to array work changed, work now has 4 dimensions instead of 1 (transpose)
array diss allocated only if required (init_3d_model)

pressure boundary condition "Neumann+inhomo" removed from the code
(check_parameters, header, poisfft, poisfft_hybrid, pres)

Errors:


bugfix: dependency added for cuda_fft_interfaces (Makefile)
bugfix: CUDA fft plans adjusted for domain decomposition (before they always
used total domain) (fft_xy)

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