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

Last change on this file since 1842 was 1834, checked in by raasch, 8 years ago

last commit documented

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