Home

Context Navigation

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |