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

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

subdomains must have identical size, i.e. grid_matching = "match" not allowed any more
parameter grid_matching removed
some obsolete variables removed

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