Home

Context Navigation

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

Last change on this file since 2372 was 2372, checked in by sward, 7 years ago
y_shift for periodic boundary conditions
Property svn:keywords set to `Id`
File size: 48.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |