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

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

last commit documented

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