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

Last change on this file since 1159 was 1159, checked in by fricke, 8 years ago

Bugfix: In case of non-cyclic lateral boundary conditions, Neumann boundary conditions for the velocity components at the outflow are in fact radiation boundary conditions using the maximum phase velocity that ensures numerical stability (CFL-condition).
Logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.

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