Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |