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

Last change on this file since 779 was 779, checked in by fricke, 12 years ago

last commit documented

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