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

Last change on this file since 1036 was 1036, checked in by raasch, 11 years ago

code has been put under the GNU General Public License (v3)

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