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

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

bugfix for mg-solver

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