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

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

2d virtual processor topology is used by default for all machines

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