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

Last change on this file since 1099 was 1093, checked in by raasch, 12 years ago

last commit documented

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