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

Last change on this file since 1086 was 1057, checked in by raasch, 12 years ago

last commit documented

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