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

Last change on this file since 4281 was 4264, checked in by scharf, 19 months ago

corrected error message string for incompatible combinations of nx, ny or nz with pdims(1) or pdims(2)

  • Property svn:keywords set to Id
File size: 50.6 KB
Line 
1!> @file init_pegrid.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: init_pegrid.f90 4264 2019-10-15 16:00:23Z schwenkel $
27! corrected error message string
28!
29! 4241 2019-09-27 06:32:47Z raasch
30! Check added to ensure that subdomain grid has at least the size as given by the number
31! of ghost points
32!
33! 4182 2019-08-22 15:20:23Z scharf
34! Corrected "Former revisions" section
35!
36! 4045 2019-06-21 10:58:47Z raasch
37! bugfix: kind attribute added to nint function to allow for large integers which may appear in
38! case of default recycling width and small grid spacings
39!
40! 3999 2019-05-23 16:09:37Z suehring
41! Spend 3 ghost points also in case of pw-scheme when nesting is applied
42!
43! 3897 2019-04-15 11:51:14Z suehring
44! Minor revision of multigrid check; give warning instead of an abort.
45!
46! 3890 2019-04-12 15:59:20Z suehring
47! Check if grid coarsening is possible on subdomain, in order to avoid that
48! multigrid approach effectively reduces to a Gauss-Seidel scheme.
49!
50! 3885 2019-04-11 11:29:34Z kanani
51! Changes related to global restructuring of location messages and introduction
52! of additional debug messages
53!
54! 3884 2019-04-10 13:31:55Z Giersch
55! id_recycling is only calculated in case of tubulent inflow
56!
57! 3761 2019-02-25 15:31:42Z raasch
58! unused variable removed
59!
60! 3655 2019-01-07 16:51:22Z knoop
61! variables documented
62!
63! Revision 1.1  1997/07/24 11:15:09  raasch
64! Initial revision
65!
66!
67! Description:
68! ------------
69!> Determination of the virtual processor topology (if not prescribed by the
70!> user)and computation of the grid point number and array bounds of the local
71!> domains.
72!> @todo: remove MPI-data types for 2D exchange on coarse multigrid level (not
73!>        used any more)
74!------------------------------------------------------------------------------!
75 SUBROUTINE init_pegrid
76 
77
78    USE control_parameters,                                                    &
79        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
80               bc_lr, bc_ns, bc_radiation_l, bc_radiation_n, bc_radiation_r,   &
81               bc_radiation_s, coupling_mode, coupling_topology, gathered_size,&
82               grid_level, grid_level_count, maximum_grid_level,               &
83               message_string, mg_switch_to_pe0_level, momentum_advec,         &
84               psolver, outflow_source_plane, recycling_width, scalar_advec,   &
85               subdomain_size, turbulent_inflow, turbulent_outflow, y_shift
86
87    USE grid_variables,                                                        &
88        ONLY:  dx
89       
90    USE indices,                                                               &
91        ONLY:  mg_loc_ind, nbgp, nnx, nny, nnz, nx, nx_a, nx_o, nxl, nxl_mg,   &
92               nxlu, nxr, nxr_mg, ny, ny_a, ny_o, nyn, nyn_mg, nys, nys_mg,    &
93               nysv, nz, nzb, nzt, nzt_mg, wall_flags_1, wall_flags_2,         &
94               wall_flags_3, wall_flags_4, wall_flags_5, wall_flags_6,         &
95               wall_flags_7, wall_flags_8, wall_flags_9, wall_flags_10
96
97    USE kinds
98     
99    USE pegrid
100   
101    USE pmc_interface,                                                         &
102        ONLY:  nested_run
103     
104    USE spectra_mod,                                                           &
105        ONLY:  calculate_spectra, dt_dosp
106
107    USE synthetic_turbulence_generator_mod,                                    &
108        ONLY:  id_stg_left, id_stg_north, id_stg_right, id_stg_south,          &
109               use_syn_turb_gen
110
111    USE transpose_indices,                                                     &
112        ONLY:  nxl_y, nxl_yd, nxl_z, nxr_y, nxr_yd, nxr_z, nyn_x, nyn_z, nys_x,&
113               nys_z, nzb_x, nzb_y, nzb_yd, nzt_x, nzt_yd, nzt_y
114
115    USE vertical_nesting_mod,                                                  &
116        ONLY:  vnested, vnest_init_pegrid_domain, vnest_init_pegrid_rank
117
118    IMPLICIT NONE
119
120    INTEGER(iwp) ::  i                        !< running index over number of processors or number of multigrid level
121    INTEGER(iwp) ::  id_inflow_l              !< ID indicating processors located at the left inflow boundary
122    INTEGER(iwp) ::  id_outflow_l             !< local value of id_outflow
123    INTEGER(iwp) ::  id_outflow_source_l      !< local value of id_outflow_source
124    INTEGER(iwp) ::  id_recycling_l           !< ID indicating processors located at the recycling plane
125    INTEGER(iwp) ::  id_stg_left_l            !< left lateral boundary local core id in case of turbulence generator 
126    INTEGER(iwp) ::  id_stg_north_l           !< north lateral boundary local core id in case of turbulence generator 
127    INTEGER(iwp) ::  id_stg_right_l           !< right lateral boundary local core id in case of turbulence generator 
128    INTEGER(iwp) ::  id_stg_south_l           !< south lateral boundary local core id in case of turbulence generator 
129    INTEGER(iwp) ::  ind(5)                   !< array containing the subdomain bounds
130    INTEGER(iwp) ::  j                        !< running index, used for various loops
131    INTEGER(iwp) ::  k                        !< number of vertical grid points in different multigrid level
132    INTEGER(iwp) ::  maximum_grid_level_l     !< maximum number of grid level without switching to PE 0
133    INTEGER(iwp) ::  mg_levels_x              !< maximum number of grid level allowed along x-direction
134    INTEGER(iwp) ::  mg_levels_y              !< maximum number of grid level allowed along y-direction
135    INTEGER(iwp) ::  mg_levels_z              !< maximum number of grid level allowed along z-direction
136    INTEGER(iwp) ::  mg_switch_to_pe0_level_l !< maximum number of grid level with switching to PE 0
137    INTEGER(iwp) ::  nnx_y                    !< quotient of number of grid points along x-direction and number of PEs used along y-direction
138    INTEGER(iwp) ::  nny_x                    !< quotient of number of grid points along y-direction and number of PEs used along x-direction
139    INTEGER(iwp) ::  nny_z                    !< quotient of number of grid points along y-direction and number of PEs used along x-direction
140    INTEGER(iwp) ::  nnz_x                    !< quotient of number of grid points along z-direction and number of PEs used along x-direction
141    INTEGER(iwp) ::  nnz_y                    !< quotient of number of grid points along z-direction and number of PEs used along x-direction
142    INTEGER(iwp) ::  numproc_sqr              !< square root of the number of processors
143    INTEGER(iwp) ::  nxl_l                    !< lower index bound along x-direction on subdomain and different multigrid level
144    INTEGER(iwp) ::  nxr_l                    !< upper index bound along x-direction on subdomain and different multigrid level
145    INTEGER(iwp) ::  nyn_l                    !< lower index bound along y-direction on subdomain and different multigrid level
146    INTEGER(iwp) ::  nys_l                    !< upper index bound along y-direction on subdomain and different multigrid level
147    INTEGER(iwp) ::  nzb_l                    !< lower index bound along z-direction on subdomain and different multigrid level
148    INTEGER(iwp) ::  nzt_l                    !< upper index bound along z-direction on subdomain and different multigrid level
149!$  INTEGER(iwp) ::  omp_get_num_threads      !< number of OpenMP threads
150
151    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  ind_all !< dummy array containing index bounds on subdomain, used for gathering
152    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxlf    !< lower index bound allong x-direction for every PE
153    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nxrf    !< upper index bound allong x-direction for every PE
154    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nynf    !< lower index bound allong y-direction for every PE
155    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nysf    !< lower index bound allong y-direction for every PE
156
157    INTEGER(iwp), DIMENSION(2) ::  pdims_remote         !< number of PEs used for coupled model (only in atmospher-ocean coupling)
158    INTEGER(iwp)               ::  lcoord(2)            !< PE coordinates of left neighbor along x and y
159    INTEGER(iwp)               ::  rcoord(2)            !< PE coordinates of right neighbor along x and y
160
161!
162!-- Get the number of OpenMP threads
163    !$OMP PARALLEL
164!$  threads_per_task = omp_get_num_threads()
165    !$OMP END PARALLEL
166
167
168#if defined( __parallel )
169
170    CALL location_message( 'creating virtual PE grids + MPI derived data types', 'start' )
171
172!
173!-- Determine the processor topology or check it, if prescribed by the user
174    IF ( npex == -1  .AND.  npey == -1 )  THEN
175
176!
177!--    Automatic determination of the topology
178       numproc_sqr = SQRT( REAL( numprocs, KIND=wp ) )
179       pdims(1)    = MAX( numproc_sqr , 1 )
180       DO  WHILE ( MOD( numprocs , pdims(1) ) /= 0 )
181          pdims(1) = pdims(1) - 1
182       ENDDO
183       pdims(2) = numprocs / pdims(1)
184
185    ELSEIF ( npex /= -1  .AND.  npey /= -1 )  THEN
186
187!
188!--    Prescribed by user. Number of processors on the prescribed topology
189!--    must be equal to the number of PEs available to the job
190       IF ( ( npex * npey ) /= numprocs )  THEN
191          WRITE( message_string, * ) 'number of PEs of the prescribed ',       &
192              'topology (', npex*npey,') does not match & the number of ',     &
193              'PEs available to the job (', numprocs, ')'
194          CALL message( 'init_pegrid', 'PA0221', 1, 2, 0, 6, 0 )
195       ENDIF
196       pdims(1) = npex
197       pdims(2) = npey
198
199    ELSE
200!
201!--    If the processor topology is prescribed by the user, the number of
202!--    PEs must be given in both directions
203       message_string = 'if the processor topology is prescribed by th' //     &
204                'e user & both values of "npex" and "npey" must be given' //   &
205                ' in the &NAMELIST-parameter file'
206       CALL message( 'init_pegrid', 'PA0222', 1, 2, 0, 6, 0 )
207
208    ENDIF
209
210!
211!-- If necessary, set horizontal boundary conditions to non-cyclic
212    IF ( bc_lr /= 'cyclic' )  cyclic(1) = .FALSE.
213    IF ( bc_ns /= 'cyclic' )  cyclic(2) = .FALSE.
214
215
216!
217!-- Create the virtual processor grid
218    CALL MPI_CART_CREATE( comm_palm, ndim, pdims, cyclic, reorder, &
219                          comm2d, ierr )
220    CALL MPI_COMM_RANK( comm2d, myid, ierr )
221    WRITE (myid_char,'(''_'',I6.6)')  myid
222
223    CALL MPI_CART_COORDS( comm2d, myid, ndim, pcoord, ierr )
224    CALL MPI_CART_SHIFT( comm2d, 0, 1, pleft, pright, ierr )
225    CALL MPI_CART_SHIFT( comm2d, 1, 1, psouth, pnorth, ierr )
226!
227!-- In case of cyclic boundary conditions, a y-shift at the boundaries in
228!-- x-direction can be introduced via parameter y_shift. The shift is done
229!-- by modifying the processor grid in such a way that processors located
230!-- at the x-boundary communicate across it to processors with y-coordinate
231!-- shifted by y_shift relative to their own. This feature can not be used
232!-- in combination with an fft pressure solver. It has been implemented to
233!-- counter the effect of streak structures in case of cyclic boundary
234!-- conditions. For a description of these see Munters
235!-- (2016; dx.doi.org/10.1063/1.4941912)
236!--
237!-- Get coordinates of left and right neighbor on PE grid
238    IF ( y_shift /= 0 ) THEN
239
240       IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
241          message_string = 'y_shift /= 0 is only allowed for cyclic ' //       &
242                           'boundary conditions in both directions '
243          CALL message( 'check_parameters', 'PA0467', 1, 2, 0, 6, 0 )
244       ENDIF
245       IF ( TRIM( psolver ) /= 'multigrid' .AND.                               &
246            TRIM( psolver ) /= 'multigrid_noopt')                              &
247       THEN
248          message_string = 'y_shift /= 0 requires a multigrid pressure solver '
249          CALL message( 'check_parameters', 'PA0468', 1, 2, 0, 6, 0 )
250       ENDIF
251
252       CALL MPI_CART_COORDS( comm2d, pright, ndim, rcoord, ierr )
253       CALL MPI_CART_COORDS( comm2d, pleft, ndim, lcoord, ierr )
254
255!
256!--    If the x(y)-coordinate of the right (left) neighbor is smaller (greater)
257!--    than that of the calling process, then the calling process is located on
258!--    the right (left) boundary of the processor grid. In that case,
259!--    the y-coordinate of that neighbor is increased (decreased) by y_shift.
260!--    The rank of the process with that coordinate is then inquired and the
261!--    neighbor rank for MPI_SENDRECV, pright (pleft) is set to it.
262!--    In this way, the calling process receives a new right (left) neighbor
263!--    for all future MPI_SENDRECV calls. That neighbor has a y-coordinate
264!--    of y+(-)y_shift, where y is the original right (left) neighbor's
265!--    y-coordinate. The modulo-operation ensures that if the neighbor's
266!--    y-coordinate exceeds the grid-boundary, it will be relocated to
267!--    the opposite part of the grid cyclicly.
268       IF ( rcoord(1) < pcoord(1) ) THEN
269          rcoord(2) = MODULO( rcoord(2) + y_shift, pdims(2) )
270          CALL MPI_CART_RANK( comm2d, rcoord, pright, ierr )
271       ENDIF
272
273       IF ( lcoord(1) > pcoord(1) ) THEN
274          lcoord(2) = MODULO( lcoord(2) - y_shift, pdims(2) )
275          CALL MPI_CART_RANK( comm2d, lcoord, pleft, ierr )
276       ENDIF
277    ENDIF
278!
279!-- Vertical nesting: store four lists that identify partner ranks to exchange
280!-- data
281    IF ( vnested )  CALL vnest_init_pegrid_rank
282
283!
284!-- Determine sub-topologies for transpositions
285!-- Transposition from z to x:
286    remain_dims(1) = .TRUE.
287    remain_dims(2) = .FALSE.
288    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dx, ierr )
289    CALL MPI_COMM_RANK( comm1dx, myidx, ierr )
290!
291!-- Transposition from x to y
292    remain_dims(1) = .FALSE.
293    remain_dims(2) = .TRUE.
294    CALL MPI_CART_SUB( comm2d, remain_dims, comm1dy, ierr )
295    CALL MPI_COMM_RANK( comm1dy, myidy, ierr )
296
297
298!
299!-- Calculate array bounds along x-direction for every PE.
300    ALLOCATE( nxlf(0:pdims(1)-1), nxrf(0:pdims(1)-1), nynf(0:pdims(2)-1),      &
301              nysf(0:pdims(2)-1) )
302
303    IF ( MOD( nx+1 , pdims(1) ) /= 0 )  THEN
304       WRITE( message_string, * ) 'x-direction: gridpoint number (',nx+1,') ', &
305                                  'is not an& integral multiple of the number',&
306                                  ' of processors (',pdims(1),')'
307       CALL message( 'init_pegrid', 'PA0225', 1, 2, 0, 6, 0 )
308    ELSE
309       nnx  = ( nx + 1 ) / pdims(1)
310    ENDIF
311
312!
313!-- Left and right array bounds, number of gridpoints
314    DO  i = 0, pdims(1)-1
315       nxlf(i)   = i * nnx
316       nxrf(i)   = ( i + 1 ) * nnx - 1
317    ENDDO
318
319!
320!-- Calculate array bounds in y-direction for every PE.
321    IF ( MOD( ny+1 , pdims(2) ) /= 0 )  THEN
322       WRITE( message_string, * ) 'y-direction: gridpoint number (',ny+1,') ', &
323                                  'is not an& integral multiple of the number',&
324                                  ' of processors (',pdims(2),')'
325       CALL message( 'init_pegrid', 'PA0227', 1, 2, 0, 6, 0 )
326    ELSE
327       nny  = ( ny + 1 ) / pdims(2)
328    ENDIF
329
330!
331!-- South and north array bounds
332    DO  j = 0, pdims(2)-1
333       nysf(j)   = j * nny
334       nynf(j)   = ( j + 1 ) * nny - 1
335    ENDDO
336
337!
338!-- Local array bounds of the respective PEs
339    nxl = nxlf(pcoord(1))
340    nxr = nxrf(pcoord(1))
341    nys = nysf(pcoord(2))
342    nyn = nynf(pcoord(2))
343    nzb = 0
344    nzt = nz
345    nnz = nz
346
347!
348!-- Set switches to define if the PE is situated at the border of the virtual
349!-- processor grid
350    IF ( nxl == 0 )   left_border_pe  = .TRUE.
351    IF ( nxr == nx )  right_border_pe = .TRUE.
352    IF ( nys == 0 )   south_border_pe = .TRUE.
353    IF ( nyn == ny )  north_border_pe = .TRUE.
354
355!
356!-- Calculate array bounds and gridpoint numbers for the transposed arrays
357!-- (needed in the pressure solver)
358!-- For the transposed arrays, cyclic boundaries as well as top and bottom
359!-- boundaries are omitted, because they are obstructive to the transposition
360
361!
362!-- 1. transposition  z --> x
363!-- This transposition is not neccessary in case of a 1d-decomposition along x
364    IF ( psolver == 'poisfft'  .OR.  calculate_spectra )  THEN
365
366       IF ( pdims(2) /= 1 )  THEN
367          IF ( MOD( nz , pdims(1) ) /= 0 )  THEN
368             WRITE( message_string, * ) 'transposition z --> x:& ',              &
369                                        'nz=',nz,' is not an integral multiple ',&
370                                        'of pdims(1)=',pdims(1)
371             CALL message( 'init_pegrid', 'PA0230', 1, 2, 0, 6, 0 )
372          ENDIF
373       ENDIF
374
375       nys_x = nys
376       nyn_x = nyn
377       nny_x = nny
378       nnz_x = nz / pdims(1)
379       nzb_x = 1 + myidx * nnz_x
380       nzt_x = ( myidx + 1 ) * nnz_x
381       sendrecvcount_zx = nnx * nny * nnz_x
382
383    ENDIF
384
385
386    IF ( psolver == 'poisfft' )  THEN
387!
388!--    2. transposition  x --> y
389       IF ( MOD( nx+1 , pdims(2) ) /= 0 )  THEN
390          WRITE( message_string, * ) 'transposition x --> y:& ',              &
391                                     'nx+1=',nx+1,' is not an integral ',     &
392                                     'multiple of pdims(2)=',pdims(2)
393          CALL message( 'init_pegrid', 'PA0231', 1, 2, 0, 6, 0 )
394       ENDIF
395
396       nnz_y = nnz_x
397       nzb_y = nzb_x
398       nzt_y = nzt_x
399       nnx_y = (nx+1) / pdims(2)
400       nxl_y = myidy * nnx_y
401       nxr_y = ( myidy + 1 ) * nnx_y - 1
402       sendrecvcount_xy = nnx_y * nny_x * nnz_y
403!
404!--    3. transposition  y --> z
405!--    (ELSE:  x --> y  in case of 1D-decomposition along x)
406       nxl_z = nxl_y
407       nxr_z = nxr_y
408       nny_z = (ny+1) / pdims(1)
409       nys_z = myidx * nny_z
410       nyn_z = ( myidx + 1 ) * nny_z - 1
411       sendrecvcount_yz = nnx_y * nny_z * nnz_y
412
413       IF ( pdims(2) /= 1 )  THEN
414!
415!--       y --> z
416!--       This transposition is not neccessary in case of a 1d-decomposition
417!--       along x, except that the uptream-spline method is switched on
418          IF ( MOD( ny+1 , pdims(1) ) /= 0 )  THEN
419             WRITE( message_string, * ) 'transposition y --> z:& ',            &
420                                        'ny+1=',ny+1,' is not an integral ',   &
421                                        'multiple of pdims(1)=',pdims(1)
422             CALL message( 'init_pegrid', 'PA0232', 1, 2, 0, 6, 0 )
423          ENDIF
424
425       ELSE
426!
427!--       x --> y
428!--       This condition must be fulfilled for a 1D-decomposition along x
429          IF ( MOD( ny+1 , pdims(1) ) /= 0 )  THEN
430             WRITE( message_string, * ) 'transposition x --> y:& ',            &
431                                        'ny+1=',ny+1,' is not an integral ',   &
432                                        'multiple of pdims(1)=',pdims(1)
433             CALL message( 'init_pegrid', 'PA0233', 1, 2, 0, 6, 0 )
434          ENDIF
435
436       ENDIF
437
438    ENDIF
439
440!
441!-- Indices for direct transpositions z --> y (used for calculating spectra)
442    IF ( calculate_spectra )  THEN
443       IF ( MOD( nz, pdims(2) ) /= 0 )  THEN
444          WRITE( message_string, * ) 'direct transposition z --> y (needed ',  &
445                                     'for spectra):& nz=',nz,' is not an ',    &
446                                     'integral multiple of pdims(2)=',pdims(2)
447          CALL message( 'init_pegrid', 'PA0234', 1, 2, 0, 6, 0 )
448       ELSE
449          nxl_yd = nxl
450          nxr_yd = nxr
451          nzb_yd = 1 + myidy * ( nz / pdims(2) )
452          nzt_yd = ( myidy + 1 ) * ( nz / pdims(2) )
453          sendrecvcount_zyd = nnx * nny * ( nz / pdims(2) )
454       ENDIF
455    ENDIF
456
457    IF ( psolver == 'poisfft'  .OR.  calculate_spectra )  THEN
458!
459!--    Indices for direct transpositions y --> x
460!--    (they are only possible in case of a 1d-decomposition along x)
461       IF ( pdims(2) == 1 )  THEN
462          nny_x = nny / pdims(1)
463          nys_x = myid * nny_x
464          nyn_x = ( myid + 1 ) * nny_x - 1
465          nzb_x = 1
466          nzt_x = nz
467          sendrecvcount_xy = nnx * nny_x * nz
468       ENDIF
469
470    ENDIF
471
472    IF ( psolver == 'poisfft' )  THEN
473!
474!--    Indices for direct transpositions x --> y
475!--    (they are only possible in case of a 1d-decomposition along y)
476       IF ( pdims(1) == 1 )  THEN
477          nnx_y = nnx / pdims(2)
478          nxl_y = myid * nnx_y
479          nxr_y = ( myid + 1 ) * nnx_y - 1
480          nzb_y = 1
481          nzt_y = nz
482          sendrecvcount_xy = nnx_y * nny * nz
483       ENDIF
484
485    ENDIF
486
487!
488!-- Arrays for storing the array bounds are needed any more
489    DEALLOCATE( nxlf , nxrf , nynf , nysf )
490
491
492!
493!-- Collect index bounds from other PEs (to be written to restart file later)
494    ALLOCATE( hor_index_bounds(4,0:numprocs-1) )
495
496    IF ( myid == 0 )  THEN
497
498       hor_index_bounds(1,0) = nxl
499       hor_index_bounds(2,0) = nxr
500       hor_index_bounds(3,0) = nys
501       hor_index_bounds(4,0) = nyn
502
503!
504!--    Receive data from all other PEs
505       DO  i = 1, numprocs-1
506          CALL MPI_RECV( ibuf, 4, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
507                         ierr )
508          hor_index_bounds(:,i) = ibuf(1:4)
509       ENDDO
510
511    ELSE
512!
513!--    Send index bounds to PE0
514       ibuf(1) = nxl
515       ibuf(2) = nxr
516       ibuf(3) = nys
517       ibuf(4) = nyn
518       CALL MPI_SEND( ibuf, 4, MPI_INTEGER, 0, myid, comm2d, ierr )
519
520    ENDIF
521
522
523#if defined( __print )
524!
525!-- Control output
526    IF ( myid == 0 )  THEN
527       PRINT*, '*** processor topology ***'
528       PRINT*, ' '
529       PRINT*, 'myid   pcoord    left right  south north  idx idy   nxl: nxr',&
530               &'   nys: nyn'
531       PRINT*, '------------------------------------------------------------',&
532               &'-----------'
533       WRITE (*,1000)  0, pcoord(1), pcoord(2), pleft, pright, psouth, pnorth, &
534                       myidx, myidy, nxl, nxr, nys, nyn
5351000   FORMAT (I4,2X,'(',I3,',',I3,')',3X,I4,2X,I4,3X,I4,2X,I4,2X,I3,1X,I3, &
536               2(2X,I4,':',I4))
537
538!
539!--    Receive data from the other PEs
540       DO  i = 1,numprocs-1
541          CALL MPI_RECV( ibuf, 12, MPI_INTEGER, i, MPI_ANY_TAG, comm2d, status, &
542                         ierr )
543          WRITE (*,1000)  i, ( ibuf(j) , j = 1,12 )
544       ENDDO
545    ELSE
546
547!
548!--    Send data to PE0
549       ibuf(1) = pcoord(1); ibuf(2) = pcoord(2); ibuf(3) = pleft
550       ibuf(4) = pright; ibuf(5) = psouth; ibuf(6) = pnorth; ibuf(7) = myidx
551       ibuf(8) = myidy; ibuf(9) = nxl; ibuf(10) = nxr; ibuf(11) = nys
552       ibuf(12) = nyn
553       CALL MPI_SEND( ibuf, 12, MPI_INTEGER, 0, myid, comm2d, ierr )       
554    ENDIF
555#endif
556
557!
558!-- Determine the number of ghost point layers
559    IF ( scalar_advec   == 'ws-scheme'  .OR.                                   &
560         momentum_advec == 'ws-scheme'  .OR.  nested_run )  THEN
561       nbgp = 3
562    ELSE
563       nbgp = 1
564    ENDIF
565
566!
567!-- Check that the number of computational grid points is not smaller than the number of
568!-- ghost points.
569    IF ( nnx < nbgp )  THEN
570       WRITE( message_string, * ) 'number of subdomain grid points along x (', nnx, ') is smaller',&
571                                  'than the number of ghost points (', nbgp, ')'
572       CALL message( 'init_pegrid', 'PA0682', 1, 2, 0, 6, 0 )
573    ENDIF
574    IF ( nny < nbgp )  THEN
575       WRITE( message_string, * ) 'number of subdomain grid points along y (', nny, ') is smaller',&
576                                  'than the number of ghost points (', nbgp, ')'
577       CALL message( 'init_pegrid', 'PA0683', 1, 2, 0, 6, 0 )
578    ENDIF
579
580!
581!-- Create a new MPI derived datatype for the exchange of surface (xy) data,
582!-- which is needed for coupled atmosphere-ocean runs.
583!-- First, calculate number of grid points of an xy-plane.
584    ngp_xy  = ( nxr - nxl + 1 + 2 * nbgp ) * ( nyn - nys + 1 + 2 * nbgp )
585    CALL MPI_TYPE_VECTOR( ngp_xy, 1, nzt-nzb+2, MPI_REAL, type_xy, ierr )
586    CALL MPI_TYPE_COMMIT( type_xy, ierr )
587
588    IF ( TRIM( coupling_mode ) /= 'uncoupled' .AND. .NOT. vnested )  THEN
589   
590!
591!--    Pass the number of grid points of the atmosphere model to
592!--    the ocean model and vice versa
593       IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
594
595          nx_a = nx
596          ny_a = ny
597
598          IF ( myid == 0 )  THEN
599
600             CALL MPI_SEND( nx_a, 1, MPI_INTEGER, numprocs, 1, comm_inter,  &
601                            ierr )
602             CALL MPI_SEND( ny_a, 1, MPI_INTEGER, numprocs, 2, comm_inter,  &
603                            ierr )
604             CALL MPI_SEND( pdims, 2, MPI_INTEGER, numprocs, 3, comm_inter, &
605                            ierr )
606             CALL MPI_RECV( nx_o, 1, MPI_INTEGER, numprocs, 4, comm_inter,  &
607                            status, ierr )
608             CALL MPI_RECV( ny_o, 1, MPI_INTEGER, numprocs, 5, comm_inter,  &
609                            status, ierr )
610             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, numprocs, 6,      &
611                            comm_inter, status, ierr )
612          ENDIF
613
614          CALL MPI_BCAST( nx_o, 1, MPI_INTEGER, 0, comm2d, ierr )
615          CALL MPI_BCAST( ny_o, 1, MPI_INTEGER, 0, comm2d, ierr ) 
616          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr )
617       
618       ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
619
620          nx_o = nx
621          ny_o = ny
622
623          IF ( myid == 0 ) THEN
624
625             CALL MPI_RECV( nx_a, 1, MPI_INTEGER, 0, 1, comm_inter, status, &
626                            ierr )
627             CALL MPI_RECV( ny_a, 1, MPI_INTEGER, 0, 2, comm_inter, status, &
628                            ierr )
629             CALL MPI_RECV( pdims_remote, 2, MPI_INTEGER, 0, 3, comm_inter, &
630                            status, ierr )
631             CALL MPI_SEND( nx_o, 1, MPI_INTEGER, 0, 4, comm_inter, ierr )
632             CALL MPI_SEND( ny_o, 1, MPI_INTEGER, 0, 5, comm_inter, ierr )
633             CALL MPI_SEND( pdims, 2, MPI_INTEGER, 0, 6, comm_inter, ierr )
634          ENDIF
635
636          CALL MPI_BCAST( nx_a, 1, MPI_INTEGER, 0, comm2d, ierr)
637          CALL MPI_BCAST( ny_a, 1, MPI_INTEGER, 0, comm2d, ierr) 
638          CALL MPI_BCAST( pdims_remote, 2, MPI_INTEGER, 0, comm2d, ierr) 
639
640       ENDIF
641 
642       ngp_a = ( nx_a+1 + 2 * nbgp ) * ( ny_a+1 + 2 * nbgp )
643       ngp_o = ( nx_o+1 + 2 * nbgp ) * ( ny_o+1 + 2 * nbgp )
644
645!
646!--    Determine if the horizontal grid and the number of PEs in ocean and
647!--    atmosphere is same or not
648       IF ( nx_o == nx_a  .AND.  ny_o == ny_a  .AND.  &
649            pdims(1) == pdims_remote(1) .AND. pdims(2) == pdims_remote(2) ) &
650       THEN
651          coupling_topology = 0
652       ELSE
653          coupling_topology = 1
654       ENDIF
655
656!
657!--    Determine the target PEs for the exchange between ocean and
658!--    atmosphere (comm2d)
659       IF ( coupling_topology == 0 )  THEN
660!
661!--       In case of identical topologies, every atmosphere PE has exactly one
662!--       ocean PE counterpart and vice versa
663          IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' ) THEN
664             target_id = myid + numprocs
665          ELSE
666             target_id = myid
667          ENDIF
668
669       ELSE
670!
671!--       In case of nonequivalent topology in ocean and atmosphere only for
672!--       PE0 in ocean and PE0 in atmosphere a target_id is needed, since
673!--       data echxchange between ocean and atmosphere will be done only
674!--       between these PEs.   
675          IF ( myid == 0 )  THEN
676
677             IF ( TRIM( coupling_mode ) == 'atmosphere_to_ocean' )  THEN
678                target_id = numprocs
679             ELSE
680                target_id = 0
681             ENDIF
682
683          ENDIF
684
685       ENDIF
686
687    ENDIF
688
689!
690!-- Store partner grid point co-ordinates as lists.
691!-- Create custom MPI vector datatypes for contiguous data transfer
692    IF ( vnested )  CALL vnest_init_pegrid_domain
693
694#else
695
696!
697!-- Array bounds when running on a single PE (respectively a non-parallel
698!-- machine)
699    nxl = 0
700    nxr = nx
701    nnx = nxr - nxl + 1
702    nys = 0
703    nyn = ny
704    nny = nyn - nys + 1
705    nzb = 0
706    nzt = nz
707    nnz = nz
708
709    ALLOCATE( hor_index_bounds(4,0:0) )
710    hor_index_bounds(1,0) = nxl
711    hor_index_bounds(2,0) = nxr
712    hor_index_bounds(3,0) = nys
713    hor_index_bounds(4,0) = nyn
714
715!
716!-- Array bounds for the pressure solver (in the parallel code, these bounds
717!-- are the ones for the transposed arrays)
718    nys_x = nys
719    nyn_x = nyn
720    nzb_x = nzb + 1
721    nzt_x = nzt
722
723    nxl_y = nxl
724    nxr_y = nxr
725    nzb_y = nzb + 1
726    nzt_y = nzt
727
728    nxl_z = nxl
729    nxr_z = nxr
730    nys_z = nys
731    nyn_z = nyn
732
733#endif
734
735!
736!-- Calculate number of grid levels necessary for the multigrid poisson solver
737!-- as well as the gridpoint indices on each level
738    IF ( psolver(1:9) == 'multigrid' )  THEN
739
740!
741!--    First calculate number of possible grid levels for the subdomains
742       mg_levels_x = 1
743       mg_levels_y = 1
744       mg_levels_z = 1
745
746       i = nnx
747       DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
748          i = i / 2
749          mg_levels_x = mg_levels_x + 1
750       ENDDO
751
752       j = nny
753       DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
754          j = j / 2
755          mg_levels_y = mg_levels_y + 1
756       ENDDO
757
758       k = nz    ! do not use nnz because it might be > nz due to transposition
759                 ! requirements
760       DO WHILE ( MOD( k, 2 ) == 0  .AND.  k /= 2 )
761          k = k / 2
762          mg_levels_z = mg_levels_z + 1
763       ENDDO
764!
765!--    The optimized MG-solver does not allow odd values for nz at the coarsest
766!--    grid level
767       IF ( TRIM( psolver ) /= 'multigrid_noopt' )  THEN
768          IF ( MOD( k, 2 ) /= 0 )  mg_levels_z = mg_levels_z - 1
769!
770!--       An odd value of nz does not work. The finest level must have an even
771!--       value.
772          IF (  mg_levels_z == 0 )  THEN
773             message_string = 'optimized multigrid method requires nz to be even'
774             CALL message( 'init_pegrid', 'PA0495', 1, 2, 0, 6, 0 )
775          ENDIF
776       ENDIF
777
778       maximum_grid_level = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
779!       
780!--    Check if subdomain sizes prevents any coarsening.
781!--    This case, the maximum number of grid levels is 1, i.e. effectively
782!--    a Gauss-Seidel scheme is applied rather than a multigrid approach.
783!--    Give a warning in this case.
784       IF ( maximum_grid_level == 1  .AND.  mg_switch_to_pe0_level == -1 )  THEN
785          message_string = 'No grid coarsening possible, multigrid ' //        &
786                           'approach effectively reduces to a Gauss-Seidel ' //&
787                           'scheme.'
788         
789          CALL message( 'poismg', 'PA0648', 0, 1, 0, 6, 0 )
790       ENDIF
791
792!
793!--    Find out, if the total domain allows more levels. These additional
794!--    levels are identically processed on all PEs.
795       IF ( numprocs > 1  .AND.  mg_switch_to_pe0_level /= -1 )  THEN
796
797          IF ( mg_levels_z > MIN( mg_levels_x, mg_levels_y ) )  THEN
798
799             mg_switch_to_pe0_level_l = maximum_grid_level
800
801             mg_levels_x = 1
802             mg_levels_y = 1
803
804             i = nx+1
805             DO WHILE ( MOD( i, 2 ) == 0  .AND.  i /= 2 )
806                i = i / 2
807                mg_levels_x = mg_levels_x + 1
808             ENDDO
809
810             j = ny+1
811             DO WHILE ( MOD( j, 2 ) == 0  .AND.  j /= 2 )
812                j = j / 2
813                mg_levels_y = mg_levels_y + 1
814             ENDDO
815
816             maximum_grid_level_l = MIN( mg_levels_x, mg_levels_y, mg_levels_z )
817
818             IF ( maximum_grid_level_l > mg_switch_to_pe0_level_l )  THEN
819                mg_switch_to_pe0_level_l = maximum_grid_level_l - &
820                                           mg_switch_to_pe0_level_l + 1
821             ELSE
822                mg_switch_to_pe0_level_l = 0
823             ENDIF
824
825          ELSE
826
827             mg_switch_to_pe0_level_l = 0
828             maximum_grid_level_l = maximum_grid_level
829
830          ENDIF
831
832!
833!--       Use switch level calculated above only if it is not pre-defined
834!--       by user
835          IF ( mg_switch_to_pe0_level == 0 )  THEN
836             IF ( mg_switch_to_pe0_level_l /= 0 )  THEN
837                mg_switch_to_pe0_level = mg_switch_to_pe0_level_l
838                maximum_grid_level     = maximum_grid_level_l
839             ENDIF
840
841          ELSE
842!
843!--          Check pre-defined value and reset to default, if neccessary
844             IF ( mg_switch_to_pe0_level < mg_switch_to_pe0_level_l  .OR.      &
845                  mg_switch_to_pe0_level >= maximum_grid_level_l )  THEN
846                message_string = 'mg_switch_to_pe0_level ' //                  &
847                                 'out of range and reset to 0'
848                CALL message( 'init_pegrid', 'PA0235', 0, 1, 0, 6, 0 )
849                mg_switch_to_pe0_level = 0
850             ELSE
851!
852!--             Use the largest number of possible levels anyway and recalculate
853!--             the switch level to this largest number of possible values
854                maximum_grid_level = maximum_grid_level_l
855
856             ENDIF
857
858          ENDIF
859
860       ENDIF
861
862       ALLOCATE( grid_level_count(maximum_grid_level),                       &
863                 nxl_mg(0:maximum_grid_level), nxr_mg(0:maximum_grid_level), &
864                 nyn_mg(0:maximum_grid_level), nys_mg(0:maximum_grid_level), &
865                 nzt_mg(0:maximum_grid_level) )
866
867       grid_level_count = 0
868!
869!--    Index zero required as dummy due to definition of arrays f2 and p2 in
870!--    recursive subroutine next_mg_level
871       nxl_mg(0) = 0; nxr_mg(0) = 0; nyn_mg(0) = 0; nys_mg(0) = 0; nzt_mg(0) = 0
872
873       nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzt_l = nzt
874
875       DO  i = maximum_grid_level, 1 , -1
876
877          IF ( i == mg_switch_to_pe0_level )  THEN
878#if defined( __parallel )
879!
880!--          Save the grid size of the subdomain at the switch level, because
881!--          it is needed in poismg.
882             ind(1) = nxl_l; ind(2) = nxr_l
883             ind(3) = nys_l; ind(4) = nyn_l
884             ind(5) = nzt_l
885             ALLOCATE( ind_all(5*numprocs), mg_loc_ind(5,0:numprocs-1) )
886             CALL MPI_ALLGATHER( ind, 5, MPI_INTEGER, ind_all, 5, &
887                                 MPI_INTEGER, comm2d, ierr )
888             DO  j = 0, numprocs-1
889                DO  k = 1, 5
890                   mg_loc_ind(k,j) = ind_all(k+j*5)
891                ENDDO
892             ENDDO
893             DEALLOCATE( ind_all )
894!
895!--          Calculate the grid size of the total domain
896             nxr_l = ( nxr_l-nxl_l+1 ) * pdims(1) - 1
897             nxl_l = 0
898             nyn_l = ( nyn_l-nys_l+1 ) * pdims(2) - 1
899             nys_l = 0
900!
901!--          The size of this gathered array must not be larger than the
902!--          array tend, which is used in the multigrid scheme as a temporary
903!--          array. Therefore the subdomain size of an PE is calculated and
904!--          the size of the gathered grid. These values are used in 
905!--          routines pres and poismg
906             subdomain_size = ( nxr - nxl + 2 * nbgp + 1 ) * &
907                              ( nyn - nys + 2 * nbgp + 1 ) * ( nzt - nzb + 2 )
908             gathered_size  = ( nxr_l - nxl_l + 3 ) * ( nyn_l - nys_l + 3 ) *  &
909                              ( nzt_l - nzb + 2 )
910
911#else
912             message_string = 'multigrid gather/scatter impossible ' //        &
913                          'in non parallel mode'
914             CALL message( 'init_pegrid', 'PA0237', 1, 2, 0, 6, 0 )
915#endif
916          ENDIF
917
918          nxl_mg(i) = nxl_l
919          nxr_mg(i) = nxr_l
920          nys_mg(i) = nys_l
921          nyn_mg(i) = nyn_l
922          nzt_mg(i) = nzt_l
923
924          nxl_l = nxl_l / 2 
925          nxr_l = nxr_l / 2
926          nys_l = nys_l / 2 
927          nyn_l = nyn_l / 2 
928          nzt_l = nzt_l / 2 
929         
930       ENDDO
931
932!
933!--    Temporary problem: Currently calculation of maxerror in routine poismg crashes
934!--    if grid data are collected on PE0 already on the finest grid level.
935!--    To be solved later.
936       IF ( maximum_grid_level == mg_switch_to_pe0_level )  THEN
937          message_string = 'grid coarsening on subdomain level cannot be performed'
938          CALL message( 'poismg', 'PA0236', 1, 2, 0, 6, 0 )
939       ENDIF
940
941    ELSE
942
943       maximum_grid_level = 0
944
945    ENDIF
946
947!
948!-- Default level 0 tells exchange_horiz that all ghost planes have to be
949!-- exchanged. grid_level is adjusted in poismg, where only one ghost plane
950!-- is required.
951    grid_level = 0
952
953#if defined( __parallel )
954!
955!-- Gridpoint number for the exchange of ghost points (y-line for 2D-arrays)
956    ngp_y  = nyn - nys + 1 + 2 * nbgp
957
958!
959!-- Define new MPI derived datatypes for the exchange of ghost points in
960!-- x- and y-direction for 2D-arrays (line)
961    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_REAL, type_x,     &
962                          ierr )
963    CALL MPI_TYPE_COMMIT( type_x, ierr )
964
965    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_REAL, type_y, ierr )
966    CALL MPI_TYPE_COMMIT( type_y, ierr )
967!
968!-- Define new MPI derived datatypes for the exchange of ghost points in
969!-- x- and y-direction for 2D-INTEGER arrays (line) - on normal grid.
970!-- Define types for 32-bit and 8-bit Integer. The 8-bit Integer are only
971!-- required on normal grid, while 32-bit Integer may be also required on
972!-- coarser grid level in case of multigrid solver.
973!
974!-- 8-bit Integer
975    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_BYTE,             &
976                          type_x_byte, ierr )
977    CALL MPI_TYPE_COMMIT( type_x_byte, ierr )
978
979    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_BYTE,                        &
980                          type_y_byte, ierr )
981    CALL MPI_TYPE_COMMIT( type_y_byte, ierr )
982!
983!-- 32-bit Integer
984    ALLOCATE( type_x_int(0:maximum_grid_level),                                &
985              type_y_int(0:maximum_grid_level) )
986             
987    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp, ngp_y, MPI_INTEGER,          &
988                          type_x_int(0), ierr )
989    CALL MPI_TYPE_COMMIT( type_x_int(0), ierr )
990
991    CALL MPI_TYPE_VECTOR( nbgp, ngp_y, ngp_y, MPI_INTEGER, type_y_int(0), ierr )
992    CALL MPI_TYPE_COMMIT( type_y_int(0), ierr )
993!
994!-- Calculate gridpoint numbers for the exchange of ghost points along x
995!-- (yz-plane for 3D-arrays) and define MPI derived data type(s) for the
996!-- exchange of ghost points in y-direction (xz-plane).
997!-- Do these calculations for the model grid and (if necessary) also
998!-- for the coarser grid levels used in the multigrid method
999    ALLOCATE ( ngp_xz(0:maximum_grid_level),                                   &
1000               ngp_xz_int(0:maximum_grid_level),                               &
1001               ngp_yz(0:maximum_grid_level),                                   &
1002               ngp_yz_int(0:maximum_grid_level),                               &
1003               type_xz(0:maximum_grid_level),                                  &
1004               type_xz_int(0:maximum_grid_level),                              &
1005               type_yz(0:maximum_grid_level),                                  &
1006               type_yz_int(0:maximum_grid_level) )
1007
1008    nxl_l = nxl; nxr_l = nxr; nys_l = nys; nyn_l = nyn; nzb_l = nzb; nzt_l = nzt
1009
1010!
1011!-- Discern between the model grid, which needs nbgp ghost points and
1012!-- grid levels for the multigrid scheme. In the latter case only one
1013!-- ghost point is necessary.
1014!-- First definition of MPI-datatypes for exchange of ghost layers on normal
1015!-- grid. The following loop is needed for data exchange in poismg.f90.
1016!
1017!-- Determine number of grid points of yz-layer for exchange
1018    ngp_yz(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1019
1020!
1021!-- Define an MPI-datatype for the exchange of left/right boundaries.
1022!-- Although data are contiguous in physical memory (which does not
1023!-- necessarily require an MPI-derived datatype), the data exchange between
1024!-- left and right PE's using the MPI-derived type is 10% faster than without.
1025    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz(0), &
1026                          MPI_REAL, type_xz(0), ierr )
1027    CALL MPI_TYPE_COMMIT( type_xz(0), ierr )
1028
1029    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz(0), ngp_yz(0), MPI_REAL, type_yz(0), &
1030                          ierr ) 
1031    CALL MPI_TYPE_COMMIT( type_yz(0), ierr )
1032
1033!
1034!-- Define data types for exchange of 3D Integer arrays.
1035    ngp_yz_int(0) = (nzt - nzb + 2) * (nyn - nys + 1 + 2 * nbgp)
1036
1037    CALL MPI_TYPE_VECTOR( nxr-nxl+1+2*nbgp, nbgp*(nzt-nzb+2), ngp_yz_int(0),   &
1038                          MPI_INTEGER, type_xz_int(0), ierr )
1039    CALL MPI_TYPE_COMMIT( type_xz_int(0), ierr )
1040
1041    CALL MPI_TYPE_VECTOR( nbgp, ngp_yz_int(0), ngp_yz_int(0), MPI_INTEGER,     &
1042                          type_yz_int(0), ierr )
1043    CALL MPI_TYPE_COMMIT( type_yz_int(0), ierr )
1044
1045!
1046!-- Definition of MPI-datatypes for multigrid method (coarser level grids)
1047    IF ( psolver(1:9) == 'multigrid' )  THEN
1048!   
1049!--    Definition of MPI-datatyoe as above, but only 1 ghost level is used
1050       DO  i = maximum_grid_level, 1 , -1
1051!
1052!--       For 3D-exchange on different multigrid level, one ghost point for
1053!--       REAL arrays, two ghost points for INTEGER arrays
1054          ngp_xz(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
1055          ngp_yz(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1056
1057          ngp_xz_int(i) = (nzt_l - nzb_l + 2) * (nxr_l - nxl_l + 3)
1058          ngp_yz_int(i) = (nzt_l - nzb_l + 2) * (nyn_l - nys_l + 3)
1059!
1060!--       MPI data type for REAL arrays, for xz-layers
1061          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz(i),       &
1062                                MPI_REAL, type_xz(i), ierr )
1063          CALL MPI_TYPE_COMMIT( type_xz(i), ierr )
1064
1065!
1066!--       MPI data type for INTEGER arrays, for xz-layers
1067          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+3, nzt_l-nzb_l+2, ngp_yz_int(i),   &
1068                                MPI_INTEGER, type_xz_int(i), ierr )
1069          CALL MPI_TYPE_COMMIT( type_xz_int(i), ierr )
1070
1071!
1072!--       MPI data type for REAL arrays, for yz-layers
1073          CALL MPI_TYPE_VECTOR( 1, ngp_yz(i), ngp_yz(i), MPI_REAL, type_yz(i), &
1074                                ierr )
1075          CALL MPI_TYPE_COMMIT( type_yz(i), ierr )
1076!
1077!--       MPI data type for INTEGER arrays, for yz-layers
1078          CALL MPI_TYPE_VECTOR( 1, ngp_yz_int(i), ngp_yz_int(i), MPI_INTEGER,  &
1079                                type_yz_int(i), ierr )
1080          CALL MPI_TYPE_COMMIT( type_yz_int(i), ierr )
1081
1082
1083!--       For 2D-exchange of INTEGER arrays on coarser grid level, where 2 ghost
1084!--       points need to be exchanged. Only required for 32-bit Integer arrays.
1085          CALL MPI_TYPE_VECTOR( nxr_l-nxl_l+5, 2, nyn_l-nys_l+5, MPI_INTEGER,  &
1086                                type_x_int(i), ierr )
1087          CALL MPI_TYPE_COMMIT( type_x_int(i), ierr )
1088
1089
1090          CALL MPI_TYPE_VECTOR( 2, nyn_l-nys_l+5, nyn_l-nys_l+5, MPI_INTEGER,  &
1091                                type_y_int(i), ierr )
1092          CALL MPI_TYPE_COMMIT( type_y_int(i), ierr )
1093
1094          nxl_l = nxl_l / 2
1095          nxr_l = nxr_l / 2
1096          nys_l = nys_l / 2
1097          nyn_l = nyn_l / 2
1098          nzt_l = nzt_l / 2
1099
1100       ENDDO
1101
1102    ENDIF
1103
1104#endif
1105
1106#if defined( __parallel )
1107!
1108!-- Setting of flags for inflow/outflow/nesting conditions.
1109    IF ( pleft == MPI_PROC_NULL )  THEN
1110       IF ( bc_lr == 'dirichlet/radiation'  .OR.  bc_lr == 'nested'  .OR.      &
1111            bc_lr == 'nesting_offline' )  THEN
1112          bc_dirichlet_l  = .TRUE.
1113       ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1114          bc_radiation_l = .TRUE.
1115       ENDIF
1116    ENDIF
1117 
1118    IF ( pright == MPI_PROC_NULL )  THEN
1119       IF ( bc_lr == 'dirichlet/radiation' )  THEN
1120          bc_radiation_r = .TRUE.
1121       ELSEIF ( bc_lr == 'radiation/dirichlet'  .OR.  bc_lr == 'nested'  .OR.  &
1122                bc_lr == 'nesting_offline' )  THEN
1123          bc_dirichlet_r  = .TRUE.
1124       ENDIF
1125    ENDIF
1126
1127    IF ( psouth == MPI_PROC_NULL )  THEN
1128       IF ( bc_ns == 'dirichlet/radiation' )  THEN
1129          bc_radiation_s = .TRUE.
1130       ELSEIF ( bc_ns == 'radiation/dirichlet'  .OR.  bc_ns == 'nested'  .OR.  &
1131                bc_ns == 'nesting_offline' )  THEN
1132          bc_dirichlet_s  = .TRUE.
1133       ENDIF
1134    ENDIF
1135
1136    IF ( pnorth == MPI_PROC_NULL )  THEN
1137       IF ( bc_ns == 'dirichlet/radiation'  .OR.  bc_ns == 'nested'  .OR.      &
1138            bc_ns == 'nesting_offline' )  THEN
1139          bc_dirichlet_n  = .TRUE.
1140       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1141          bc_radiation_n = .TRUE.
1142       ENDIF
1143    ENDIF
1144!
1145!-- In case of synthetic turbulence geneartor determine ids.
1146!-- Please note, if no forcing or nesting is applied, the generator is applied
1147!-- only at the left lateral boundary.
1148    IF ( use_syn_turb_gen )  THEN
1149       IF ( bc_dirichlet_l )  THEN
1150          id_stg_left_l = myidx
1151       ELSE
1152          id_stg_left_l = 0
1153       ENDIF
1154       IF ( bc_dirichlet_r )  THEN
1155          id_stg_right_l = myidx
1156       ELSE
1157          id_stg_right_l = 0
1158       ENDIF
1159       IF ( bc_dirichlet_s )  THEN
1160          id_stg_south_l = myidy
1161       ELSE
1162          id_stg_south_l = 0
1163       ENDIF
1164       IF ( bc_dirichlet_n )  THEN
1165          id_stg_north_l = myidy
1166       ELSE
1167          id_stg_north_l = 0
1168       ENDIF
1169
1170       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1171       CALL MPI_ALLREDUCE( id_stg_left_l, id_stg_left,   1, MPI_INTEGER,       &
1172                           MPI_SUM, comm1dx, ierr )
1173
1174       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1175       CALL MPI_ALLREDUCE( id_stg_right_l, id_stg_right, 1, MPI_INTEGER,       &
1176                           MPI_SUM, comm1dx, ierr )
1177
1178       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1179       CALL MPI_ALLREDUCE( id_stg_south_l, id_stg_south, 1, MPI_INTEGER,       &
1180                           MPI_SUM, comm1dy, ierr )
1181
1182       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1183       CALL MPI_ALLREDUCE( id_stg_north_l, id_stg_north, 1, MPI_INTEGER,       &
1184                           MPI_SUM, comm1dy, ierr )
1185
1186    ENDIF
1187 
1188!
1189!-- Broadcast the id of the inflow PE
1190    IF ( bc_dirichlet_l )  THEN
1191       id_inflow_l = myidx
1192    ELSE
1193       id_inflow_l = 0
1194    ENDIF
1195    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1196    CALL MPI_ALLREDUCE( id_inflow_l, id_inflow, 1, MPI_INTEGER, MPI_SUM, &
1197                        comm1dx, ierr )
1198
1199!
1200!-- Broadcast the id of the recycling plane
1201!-- WARNING: needs to be adjusted in case of inflows other than from left side!
1202    IF ( turbulent_inflow ) THEN
1203   
1204       IF ( NINT( recycling_width / dx, KIND=idp ) >= nxl  .AND.                                   &
1205            NINT( recycling_width / dx, KIND=idp ) <= nxr )  THEN
1206          id_recycling_l = myidx
1207       ELSE
1208          id_recycling_l = 0
1209       ENDIF
1210       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1211       CALL MPI_ALLREDUCE( id_recycling_l, id_recycling, 1, MPI_INTEGER, MPI_SUM, &
1212                           comm1dx, ierr )
1213                           
1214    ENDIF
1215
1216!
1217!-- Broadcast the id of the outflow PE and outflow-source plane
1218    IF ( turbulent_outflow )  THEN
1219
1220       IF ( bc_radiation_r )  THEN
1221          id_outflow_l = myidx
1222       ELSE
1223          id_outflow_l = 0
1224       ENDIF
1225       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1226       CALL MPI_ALLREDUCE( id_outflow_l, id_outflow, 1, MPI_INTEGER, MPI_SUM, &
1227                           comm1dx, ierr )
1228
1229       IF ( NINT( outflow_source_plane / dx ) >= nxl  .AND. &
1230            NINT( outflow_source_plane / dx ) <= nxr )  THEN
1231          id_outflow_source_l = myidx
1232       ELSE
1233          id_outflow_source_l = 0
1234       ENDIF
1235       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1236       CALL MPI_ALLREDUCE( id_outflow_source_l, id_outflow_source, 1, &
1237                           MPI_INTEGER, MPI_SUM, comm1dx, ierr )
1238
1239    ENDIF
1240
1241    CALL location_message( 'creating virtual PE grids + MPI derived data types', 'finished' )
1242
1243#else
1244    IF ( bc_lr == 'dirichlet/radiation' )  THEN
1245       bc_dirichlet_l = .TRUE.
1246       bc_radiation_r = .TRUE.
1247    ELSEIF ( bc_lr == 'radiation/dirichlet' )  THEN
1248       bc_radiation_l = .TRUE.
1249       bc_dirichlet_r = .TRUE.
1250    ENDIF
1251
1252    IF ( bc_ns == 'dirichlet/radiation' )  THEN
1253       bc_dirichlet_n = .TRUE.
1254       bc_radiation_s = .TRUE.
1255    ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
1256       bc_radiation_n = .TRUE.
1257       bc_dirichlet_s = .TRUE.
1258    ENDIF
1259#endif
1260
1261!
1262!-- At the inflow or outflow, u or v, respectively, have to be calculated for
1263!-- one more grid point.
1264    IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
1265       nxlu = nxl + 1
1266    ELSE
1267       nxlu = nxl
1268    ENDIF
1269    IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
1270       nysv = nys + 1
1271    ELSE
1272       nysv = nys
1273    ENDIF
1274
1275!
1276!-- Allocate wall flag arrays used in the multigrid solver
1277    IF ( psolver(1:9) == 'multigrid' )  THEN
1278
1279       DO  i = maximum_grid_level, 1, -1
1280
1281           SELECT CASE ( i )
1282
1283              CASE ( 1 )
1284                 ALLOCATE( wall_flags_1(nzb:nzt_mg(i)+1,         &
1285                                        nys_mg(i)-1:nyn_mg(i)+1, &
1286                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1287
1288              CASE ( 2 )
1289                 ALLOCATE( wall_flags_2(nzb:nzt_mg(i)+1,         &
1290                                        nys_mg(i)-1:nyn_mg(i)+1, &
1291                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1292
1293              CASE ( 3 )
1294                 ALLOCATE( wall_flags_3(nzb:nzt_mg(i)+1,         &
1295                                        nys_mg(i)-1:nyn_mg(i)+1, &
1296                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1297
1298              CASE ( 4 )
1299                 ALLOCATE( wall_flags_4(nzb:nzt_mg(i)+1,         &
1300                                        nys_mg(i)-1:nyn_mg(i)+1, &
1301                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1302
1303              CASE ( 5 )
1304                 ALLOCATE( wall_flags_5(nzb:nzt_mg(i)+1,         &
1305                                        nys_mg(i)-1:nyn_mg(i)+1, &
1306                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1307
1308              CASE ( 6 )
1309                 ALLOCATE( wall_flags_6(nzb:nzt_mg(i)+1,         &
1310                                        nys_mg(i)-1:nyn_mg(i)+1, &
1311                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1312
1313              CASE ( 7 )
1314                 ALLOCATE( wall_flags_7(nzb:nzt_mg(i)+1,         &
1315                                        nys_mg(i)-1:nyn_mg(i)+1, &
1316                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1317
1318              CASE ( 8 )
1319                 ALLOCATE( wall_flags_8(nzb:nzt_mg(i)+1,         &
1320                                        nys_mg(i)-1:nyn_mg(i)+1, &
1321                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1322
1323              CASE ( 9 )
1324                 ALLOCATE( wall_flags_9(nzb:nzt_mg(i)+1,         &
1325                                        nys_mg(i)-1:nyn_mg(i)+1, &
1326                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1327
1328              CASE ( 10 )
1329                 ALLOCATE( wall_flags_10(nzb:nzt_mg(i)+1,        &
1330                                        nys_mg(i)-1:nyn_mg(i)+1, &
1331                                        nxl_mg(i)-1:nxr_mg(i)+1) )
1332
1333              CASE DEFAULT
1334                 message_string = 'more than 10 multigrid levels'
1335                 CALL message( 'init_pegrid', 'PA0238', 1, 2, 0, 6, 0 )
1336
1337          END SELECT
1338
1339       ENDDO
1340
1341    ENDIF
1342
1343 END SUBROUTINE init_pegrid
Note: See TracBrowser for help on using the repository browser.