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

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

last commit documented

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